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ABSTRACT 

A failure detection filter is applied to the detec- 
tion of actuator and sensor failures on a free-free beam. 
Computer simulation tests are used to verify the filter 
design and study the effect of unmodeled modes on filter 
performance . 

In actuator tests, the failure signal to spillover 
noise ratio was found to be greatest when the filter 
bandwidth was 5 rad/sec beyond the input frequency. Ob- 
servation spillover, however, was found to vary widely in 
tests run under similar conditions (same input frequency 
and filter poles) but with different detector gains. 

In sensor tests, the maximum signal-to-noise ratio 
for varying filter bandwidth depended upon the initial 
conditions placed on the unmodeled modes; the performance 
was good even for initial amplitudes on the first un- 
modeled mode 7.5% of that on the last modeled mode. 

Data-sampling tests on filters designed for continu- 
ous data processing but employed in a sampled data mode 
revealed that adequate filter performance could be 
achieved only when the sampling rate was considerably 
beyond the natural frequency of the last system mode. 
Stability problems were encountered when the filter band- 
width became too high relative to the sampling rate. 

The failure simulation tests suggest high sampling 
rates and sensor post-filtering to deal with the problems 
posed by sampling phase lag and observation spillover. 
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CHAPTER I 


INTRODUCTION 

With the advent of the space shuttle, aerospace engineers 
are contemplating the assembly and deployment of some very 
large space structures. Some structures under consideration ' 
include antennas and reflectors 100 meters in diameter and 
solar power satellites as large as 20 x 10 kilometers. Un- 
like the spacecraft of previous decades, these large struc- 
tures have little inherent rigidity due to their low mass and 
large size. If the natural damping is not somehow increased, 
periodic disturbances such as gravity gradient and solar pres- 
sure which are close to the low natural frequencies of the 
structure will cause large dynamic overstresses that will 
eventually tear the structure apart. 

The solar power satellite provides a good example of the 
types of overwhelming issues one would typically encounter in 
designing a control system for a large space structure. In 
order to adequately damp the many vibrational modes of the 
satellite, hundreds of thrusters and control moment gyros may 
be required to supplement passive damping. The system designer 
will have to decide how many actuators and sensors to use and 
where to place them on the structure. For example, rate gyro 
sensors and control moment gyros could be located almost any- 
where on a truss-like structure. The control engineer will 
then have to decide what kind of control law to implement in 
order to maintain satisfactory structural rigidity. Obviously, 


7 



the control system cannot incorporate all the structural 
modes in its model, so care must be taken when controlling 
the disturbance-induced vibrations in the low frequency 
modes that the control does not spillover into the higher 
frequency unmodeled modes. 

One factor which should not be overlooked in either 
the design or operation of the control system is the likeli- 
hood of some failures among actuators and sensors. For ex- 
ample, if the interval between maintenance visits is three 
years and the control system utilizes a total of 400 sensors 
and actuators each with an exponential distribution of time 
to failure with a mean time to failure of 100,000 hours, the 
expected number of failures in this interval is 92, and the 
probability that there will be no failures is 2 x 10~^^. 

Thus even with a very optimistic mean time to failure, it 
is virtually certain that failures will occur. 

One of the major issues in dealing with component un- 
reliability in control systems is how to detect a failure 
and identify the failed component. This thesis is concerned 
with one method of doing failure detection and identification 
(FDD . 

Many approaches to FDI have been used, the simplest of 
which involves triplication of components; a discrepancy be- 
tween the signals of two like sensors signifies a failure, 
and comparison with the third determines which of the two has 
failed. Though simple, this method rapidly becomes costly and 
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even bulky for certain applications. 

There are several approaches to FDI which require 
specification of failure inodes ahead of time, but one which 
does not is generalized parity relations.^ This method uses 
sensor data from several time steps to detect failures 
rather than data from duplicate sensors at the same time 
instant. This approach has the obvious advantage of re- 
quiring fewer components, but it turns out to be very 
susceptible to plant disturbances and sensor noise. This 
detection routine also performs poorly when there is model 
error present, whether it be in the form of modal trunca- 
tion or frequency errors. 

A closed-loop method, the failure detection filter, 
can simultaneously monitor many different types of components, 
including sensors, actuators, and dynamic elements of the 
system. As with any other observer, the detection filter 
incorporates a linear -dynamic model of the system to esti- 
mate the true states of the system. Since the model re- 
ceives the same control inputs as the true system, the out- 
puts of the system and filter will normally match resulting 
in an output error of zero. However, when a component fails, 
the output error will no longer be zero, signifying that a 
failure has occurred. The failed component can then be 
identified by the fixed line or plane to which the output 
error is restricted by the detection filter. 

The failure detection filter was first proposed by 
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Beard (1) in 1971 for deterministic systems. The theory 
was later expanded by Jones (2) to stochastic and sampled- 
data systems. Though not strictly valid for sampled-data 
systems, the detection filter will behave satisfactorily 
for sufficiently high sampling rates. Besides the applica- 
tion by Jones to a lateral mode autopilot, the failure de- 
tection filter was applied by VanderVelde (8) and Gerard 
(9) to the computer control of a guideway vehicle, and by 
Meserole (3) to fault-tolerant control of a turbofan engine. 
In neither of these previous applications was the filter 
designed to detect a sensor failure when the sensor output 
was not measuring a state directly. This is also the first 
time model error has been introduced into the filter. 

The next chapter summarizes the main concepts of 
failure detection theory along with an analytic design pro- 
cedure for the filter. Chapter 3 proposes a computational 
design procedure for the deterministic filter based mainly 
upon algorithms siaggested in Appendix A of Beard (1) . In 
Chapter 4, simulation results of actuator and sensor failures 
in deterministic systems are presented along with some re- 
sults on data-sampling. Finally, some conclusions are out- 
lined in Chapter 5. 


V 


10 



CHAPTER II 


FAILURE DETECTION FILTER THEORY 

Detection filter theory is based upon vector-space con- 
cepts involving the state estimation errors generated by the 
filter following component failures. The major feature of 
the failure detection filter is that the output error is small 
while the system is functioning normally, and following the 
failure of a system component that error is significantly 
larger and appears only in a single direction or plane — ^that 
direction or plane indicating which component has failed. Thus 
the filter provides the basis for both detection of component 
failures and isolation of the faulty component. It is not 
necessary to specify in advance the possible modes of component 
failures . 

In this chapter, the structure of the failure detection 
filter is first presented along with failure models for both 
actuator and sensor malfunctions and plant dynamics changes. 

The concept of failure "detectability" will then be introduced 
followed by the filter design theory for both fully measurable 
systems (rank C=n) and partially measurable systems (rank C<n) . 
Finally, a two-mode beam example is used to illustrate the 
analytic design procedure. 

2.1 Detection Filter Structure 

The structure of the failure detection filter depicted 
in Fig. 2.1 is the same as that of any other state estimation 
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Fig. 2.1; Failure detection filter block diagram. 

(taken from reference [1] ) 
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filter monitoring a linear time invariant system. The differ- 
ence is in the feedback gain D which is designed so that the 
output error will be held to a fixed direction or plane depend- 
ing upon what component has failed. If D were the Kalman gain, 
the filter would minimize the mean square error between the 
states x(t) and the filter estimates x(t) when noise is intro- 
duced into the system. For the case of an observer, $(t) ap>- 
proaches x(t) asymptotically since no noise is present. 

The system being monitored by the detection filter must 
be linear, time-invariant, and observable. The consequences 
of unobservability will be pointed out as they are encountered. 
In this work, mainly continuous inputs and measurements are 
considered. 

The system is represented by the linear equations: 

Axio+'BuW 

= t xUl 

where the state vector x(t) has dimension n, the control vector 
u(t) has m inputs, and the measurement vector y(t) has p outputs. 

The filter incorporates the model: 

AxW t'BuU') ( 2 . 1 - 2 ) 


where x is the filter state vector and y^(t) the predicted 
measurement vector. 

The state and measurement error dynamics are 


4W= [A-T)Cl eW 

CeW 


13 



where e (t) - x (t) - x(t) and £ (t) = ^(t) - ^(t) . If the 
eigenvalues of A-DC have negative real parts, the steady 
state error will be zero so x(t) will track x(t) . Equation 
2.1-3 is valid, however, only when ^A,B,c”^ of the filter 
in (2.1-2) match that of the true system represented by 
(2.1-1) . 

2.2 Failure Models 


As mentioned above, (2.1-3) is valid only when the trip- 
lets [a,B,c| characterizing the system and filter match. If 
there is a component failure or change in plant dynamics, 
(2.1-1) will no longer represent the true system. 

Consider an actuator j that has either failed completely 
or gone awry. The control vector, u(t) , is given by 

U.CO = (2.2- 

A 

where u (t) is the commanded input (and filter input) , e is 

th ^ 

a unit vector in the j direction of dimension m, and n(t) 

is a scalar time' function depending upon the type of failure. 

The new state equation becomes 


X CO = Ax CO u^CO + ti nCO 


( 2 . 2 - 2 ) 


where bj is the column of B corresponding to the j actuator 
and f\(t) is some scalar time function. If the actuator fails 
completely in the off mode, r\(t) is simply minus the in- 

put, -u (t) , and the control effectiveness matrix for the true 

^ th 

system is the same as for the filter except that the j column 


of B is deleted. 
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Subtracting (2.1-2) from (2.2-2), the error equations be- 
come 

lA-DCTleW+tinW (2.2-3) 

tW' CeCV) 

The vector b^ is called the "event vector" because it is the 
driving vector in the event of failure of the actuator. 

Detection filter theory will show that for any event vector 
it is possible to find a D such that €(t) maintains a 
fixed direction in response to fn(t). 

Similarly, consider a sensor failure represented by the 
new measurement relation 

Cx Ct') + (2.2-4) 

• th 

where e is a unit vector in the j direction and n(t) is 

an arbitrary scalar function of time. Complete failure of 

the sensor in the zero-output mode is modeled by letting 

T ' 

n(t) = ~ Cj x(t) , the predicted sensor output, and the measure- 
ment matrix for the true system is that of the filter with the 
row deleted. 

The new error equations are generated analogously with 
(2.2-4) replacing ^(t) in (2.1-1); 

eCO= [VDC.1 (2 2-5) 

Unlike the actuator case, the output error in the sensor failure 
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mode will have two directions; e due directly to the sensor 

failure and another component caused by the measurement error 

effect in the filter feedback loop. Even though D can be 

chosen to make t(t) due to dj unidirectional, this direction 

A 

will usually not be e , and therefore the output error will 

span a plane rather than lie in a fixed direction. 

Finally, let us consider changes in plant dynamics. Let 
th. 

the ij element of A, for instance, change by the amount 
Aa^j . The new state equation is ; 

-v'BaCl'i -v e«j (2.2-6) 

This model represents the effect of an alteration in the deriva- 
tive of x^(t) due to dynamics involving Xj (t) . The error equa- 
tions become 

ew- [^-Dcle^O+ AaiiXjtoL; (2-2-7) 

C.%w 

Comparing this, with Eq. (2.2-3), we see that the two pairs of 
equations have the same form with the scalar time function in 
this case being Aa^^x^ (t) and the new event vector as the 
n-dimensional i unit vector. 

Actuator and plant dynamics failure models are, in fact, 
called "input failure models" since both can be represented as 
extraneous additive inputs to the system. tVhen the failure 
term also appears in the measurement equation [as in Eq. (2.2-5) 
for sensor failures] , the error equations are referred to as 
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the "measurement failure models." Only the actuator input 
and measurement failures will be considered in this work. 

2.3 Detection Filter Design 

The preceding discussion of failure models motivates the 
following definition of "detectability" found in Beard [1] : 
The event associated with the vector ^ in the state error 
equation 

rK-t)C.l%W + InK) , 2 . 3-1 

is "detectable" if there exists a gain matrix D such that; 

maintains a fixed direction in the output space 
(ep(t) is the settled-out solution of Eq. 2.3-1) 

(ii) All the eigenvalues of (A-DC) can be arbitrarily speci- 
fied (this condition may be relaxed for fixed but stable 
eigenvalues) . 

If ^ is detectable, then one can identify the failed com- 
ponent by checking the direction of the output error. For a 
measurement failure, the definition of detectability must be 
revised (see Section 2.3.3) and the output error is confined 
to a plane. 

In the next several sections, the design of D for both 
fully measurable systems (rank C=n) and partially measurable 
systems (rank C^) will be presented. Finding D for the first 
case is simple since the state vector x(t) can be solved 
uniquely given y^(t) . However, when the state vector is not 
fully measurable, the design of D is much more difficult and 
generates the need for a more advanced methodology. 
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2.3.1 Fully Measurable Systems 

Let NS be the mamber of sensors or rows of C. To satis- 
fy condition (ii) of detectability, choose A-DC to equal - Cl 
where C is a positive scalar constant. Then the eigenvalues 
of A-DC are all - c and D is given uniquely by 

"D= (.K+ <rl') C"' (2.3-2) 

for rankC=n=NS and nonuniquely as 

"D’ CK-vtrIHcrC.y' C’' (2.3-3) 


for rankC=n<NS. 

For an actuator failure, £ = b and the solution of Eg. 
(2.2-3) is 


eu') - Ji ecto') 


(2.3-4) 


The first term on the right is due to the initial conditions 
of the filter and will asymptotically approach zero for CJ* > 0. 
The second term is the driven or settled solution — ^it comes 
about because the filter models the i^^ actuator as working 
properly (u(t) = u (t) ) whereas in actuality it is not. 

Since b. is constant, the settled-out solution is 
—1 

* “b ; f r\{X) (2.3-5) 


Therefore, e (t) lies in the direction b. and a settled output 

“*S X, 

error in the direction Cb^^ indicates a malfunction of the i^^ 
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actuator. 


Now consider a failure in the i^^ sensor as modeled by 
(2.2-4). Retaining the same choice for A-DC, the solution 
to (2.2-5) is 











(2.3-6) 


The settled output error is given by 


to 


/s. 

Therefore, an output error lying in the plane spanned by CDe . 

A th 

and e . indicates a failure of the i sensor. 

-pi 

Thus we see that finding D in the fully measurable case 
is indeed trivial and given by (2.3-2) for A-DC = -0*1. Choos- 
ing A-DC = - O' I enabled the assignment of eigenvalues directly 
and also produced readily identifiable output error directions 
for component malfunctions because of the simple form of the 
state error transition matrix. 

2.3.2 Partially Measurable Systems 

Because C is not invertible in partially measurable sys- 
tems, D cannot be solved by a simple relation like (2.3-2). What 
this means physically is that when an event occurs, an actuator 
failure for instance, the output error which results will not 
be uniquely associated with the state error caused by the 
event. Thus since the output error is the only accessible 
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signal, D will not be unique in satisfying the conditions of 
detectability. It will therefore be necessary to explore 
the relation of D to the space spanned by state errors during 
an event. 


2.3. 2.1 Detection Generator 


Consider the system equations (2.1-1) for the case of a 
single actuator: 


+-felUW (. 2 . 2 - 1 ) 

From linear systems theory, the controllable space of u(t) , 

{the space of states which are possible to reach from the 
origin using u(t)] is the range space of 


U,i • • • : K”-'y 


(2.3-8) 


The event vector ^ has a similar meaning with respect to 
the error space e(t) of (2.3-1). In this case, the error 
space that can be reached in the event of a failure is the 
range space of 

' “[{iOV-DCHi ■ ■ ■ i (2.3-9) 

or the "detection space" of The output error space there- 

fore spans CW^, and condition (i) of detectability is satisfied 
only if the rank of CW^ is 1. Any matrix D which satisfies 
this condition is called a "detector gain" for £. 

Associated with f is another vector also in defined 
by the following theorem proven in Beard: 
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Detection generator theorem - 


If (i) (A,C) is an observable pair 

(ii) rank = k 
(iii) rank CW^ = 1 

Then there exists an n-vector g. such that 


and 


L t i CM 


g ^ 0 



(2.3-10) 


The latter relation in (2.3-10) guarantees that the controllable 
space of 2^ with respect to (A-DC) is of dimension k and so 
matches the controllable space of Note also that the first 

equation in (2.3-10) yields 


[ (K-DCVl-- •UK-'DC.f'g] 
= t %’l • i 


(2.3-11) 


SO the set of vectors Ag^, ...,A^ ^g| form a basis for the 

controllable space of f . 

Since £ is the generator of W^, it can be expressed as 


f - . .-V o(y^ 


(2.3-12) 


where the are scalars. Premultiplying by C and using 
(2.3-10) , 

Cf = o(vc (2.3-13) 

If C£ 0, then f 0, and since the magnitude of g_ is not 
restricted by (2.3-10), we are free to choose = 1. 

If in general, for some nonnegative integer 
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C KH ' q 

tK*'i 0 

then 

O 

and £ is chosen so that The fact that (A,C) is an 

observable pair will guarantee that (2.3-14) holds for 
juii: k-1. Using (2.3-15) in (2.3-12) yields 

. f= + . . .+ (2.3-16) 

where ^ satisfies (2.3-14). 

The n-vector g which satisfies (2.3-10) and (2.3-16) is 
called a order "detection generator" for The detection 

generator theorem has shown that as long as (A,C) is an ob- 
servable pair, there always exists a detection generator as- 
sociated with a detector gain. 

2. 3. 2. 2 Detector Gain 

The preceding development of the detection generator serves 
as a prelude to the following theorem from which D can be found: 
Detector gain theorem — 

If (i) (A,C) is an observable pair 

(ii) rank = k 
(iii) rank CW^ = 1 

and k eigenvalues of [A-DC] associated with 
are given by the roots of 




(2.3-14) 


for i = O, . . . , )JL-\ 


(2.3-15) 
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(2.3-17) 


5 ^ -V + . . . -V- -p, -- o 

Then D is a solution of 

DC./^ 'i - 'Pv^ t g V (2.3-18) 

where are scalars and is a order detection 

generator for ^ 

Conversely, any D satisfying (2.3-18) is a detector gain 
for £ and k eigenvalues of [A-DC] associated with 
are given by the roots of (2.3-17). 

This theorem is perhaps the main result of detection filter 
theory and follows mainly from the Cayley-Hamilton Theorem applied 
to [A-DC] and invoking (2.3-11). 

Premultiplying (2.3-16) by CA^ and using (2.3-10) yields 

^ (2.3-19) 

which when substituted into (2.3-18) gives 

( 2 . 3 - 20 ) 

where JJ.= 0 corresponds to the condition C^ ^ 0. The rest of 
the results will be derived for 0, but the corresponding 

general results can be obtained by simply replacing ^ by A^^. 

The general solution of (2.3-20) is expressible in the 

form 

T) = 

where (2.3-21) 

Di,= iciV 

Dh - T)' li- tiKtiycy' iciY] 

with 2^ being the right side of (2.3-20) and D' an arbitrary 
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matrix the same dimensions as D. Dp is the particular solution 
of D and satisfies (2.3-20). Therefore, Dp is alone sufficient 
to guarantee condition (i) of detectability and specify k 
eigenvalues of [A-DC] . Djj is the homogeneous solution of D, 
i.e., D„Cf = 0, and represents the freedom left in D after hav- 
ing satisfied (2.3-20). Condition (ii) of detectability will 
be satisfied only if the remaining n-k eigenvalues of [A-DC] 
can be assigned using D*. 

2. 3. 2. 3 Detection Space 


Beard proposes the following auxiliary lemma to show that 
the remaining eigenvalues of [A-DC] can be assigned (the actual 
assignment of these remaining eigenvalues is discussed in Chapter 
3) : 


Auxiliary detector gain lemma — 

For A' , C' , and D' given by 

= , 2 . 3 . 22 , 

C' = C - CVUctVCiT' (tiVC 

D' 

The number of eigenvalues of [A'-D'C'] = [A-DC] that can 
be specified by D' is q' = rank M' where 

ic i CK'i •• -l C.'K"“r (2.3-23) 

This lemma shows that the total number of eigenvalues of [A-DC] 
that can be specified while constraining D to be a detector gain 
is k + q'. Therefore, condition (ii) of detectability will be 
satisfied only if n = k + q'. 

The next lemma shows that M' does not depend upon the par- 
ticular detection generator used or its order k even though it 
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appears to through A' : 


Detection space theorem — 

For K defined by 

K’ K - KiUC.tyC.iV (C£rc 

M' is also equivalent to 

M\'= [ t'\ C'Kl • • • i CK"'']"" (2.3-24) 

since C and K depend only on A, C, and M' is the same 
regardless- of what detection generator is used. Thus the number 
of eigenvalues of [A-DC] that can be specified while constraining 
D to be a detector gain is always the same no matter what g. is 
used to find Dp. The null space of M' is called the "detection 
space" of f and corresponds to the controllable space of ^ with 
respect to [A-DC] . The dimension of the detection space is 
called the "detection order" V , and a detection generator of 
order V is called the "maximal generator" for 

The rank of the null space of M' is V= n - q', so if q' 
eigenvalues of [A-DC] are assignable using D' , then a detection 
generator of order k = V is needed to assign all the eigenvalues 
and satisfy condition (ii) of detectability. Beard shows that if 
(A,C) is an observable pair, £ will have a maximal generator and 
it is unique. Thus, the event ^ is guaranteed detectable if 
(A,C) is observable. Beard also shows that although observabil- 
ity is necessary to satisfy condition (ii) of detectability, it 
is not necessary to satisfy condition (i) if ^ is in the observable 
space. The number of assignable eigenvalues in this case is the 
rank of the observability matrix M defined by (2.3-23) without 
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the primes. Therefore, if one is willing to accept n-rank M 
unassignable eigenvalues, every observable event is detectable 
in this weaker sense for the unobservable pair (A,C) . 

The following theorem summarizes the conditions for detecta- 
bility : 

Detectability Theorem — 

Every vector in the state space is detectable if and 
only if (A,C) is an observable pair. 

2.3.3 Sensor Detectability 

In deriving (2.2-5), we noted a fundamental difference be- 
tween measurement and input failures; whereas an actuator failure 
only drives the state error through the event vector b^ , a sensor 
failure drives the state error through the output feedback term 
as well as the event vector d ^ . Consequently, the best the filter 
can do is constrain the output error to a two-dimensional plane. 
Even though this was shown only for the fully measurable case. 

Beard proves this is true for partially measurable systems as well. 

Because the event vector d^ is at the discretion of the de- 
signer and is not related at all to the sensor itself, the theorems 
in the previous sections are not directly applicable. Instead, 
the i^^ row of C will be associated with a sensor failure just as 
the i^^ column of B was associated with an actuator failure. This 
motivates the following modified definition of detectability for 
sensor failures; 

j.V| ~J~ 

The i row of C, C; , is "sensor- detectable" if there 
exists a gain matrix D such that; 
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(i) §p(t) is constrained to a 2D-plane (where is 

the settled-out solution of (2.2-5)) 

(ii) the eigenvalues of [A-DC] can be specified arbitrarily 
The following theorem is the sensor counterpart to the 
Detectability Theorem: 

Sensor Detectability Theorem — 

rn 

If (A,C) is an observable pair and is linearly inde- 
pendent of all the other rows of C, then is sensor de- 

tectable. 

The filter design procedure for sensors is stated in two 
steps which also serve as a basis for the proof of this theorem: 

(1) first choose D to be a detector gain for by making 
it a solution of (2.3-20), where 

~ (2.3-25) 

(2) then make D' a detector gain for d^^ with respect to 
the pair (A* ,C' ) 

Before proving the theorem, we first outline two implica- 
tions of (2.3-25). First, the existence of ^ is guaranteed be- 

T - ■ 

cause c^ is linearly independent of the other rows of C (pseudo- 
inverse of C exists). Therefore by the detectability theorem, f 
is detectable. Second, d^ is now fixed because (2.3-20) becomes 

4i- -p,^+ . . .+-pvK'‘g -I- (2.3-26) 

where V is the detection order of £, and g. is its maximal 

generator. 

Assume for the moment that d^ lies in the observable space 
of (A',C*) and prove that condition (i) of sensor detectability 
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holds. This assumption plus the detectability theorem guarantee 

that the detector gain D' exists for the event d^ .with respect 

to (A',C). Furthermore, C'e'g(t) lies in a fixed direction 

where e' (t) is the settled-out solution of 
— s 

e W - (,K--D'C") e'W - 4; nW ( 2 . 3 - 27 ) 

Since A'-D'C is equal to A-DC, ^ (t) is the settled state error 
for the pair (A,C) using the gain D as well. This fixed direc- 
tion can therefore be expressed as 

C'e',C•t^= ’ (2.3-28) 

where y , is this fixed direction and m, (t) is a scalar function 
=i-d d 

depending upon n(t). 

From (2.3-22), C' is given by 

t'=C-|,,aT 

so that 

^ Ip; C\ IsW 

From (2.2-5), the output error is 

Ip.tnCtU C 


(2.3-29) 


(2.3-30) 


(2.3-31) 


Since the quantity in brackets is a scalar, ^(t) lies in the plane 

spanned by ^ and ^ . Beard shows that ^ can be broken up into 

0 

a component in the e direction and the CK'^Af direction, where K 


-Pi 


is given in (2.3-24) and is the lowest integer such that 
CK*^Af ^ 0. Thus the output error lies in the plane formed by 


e and CK'^Af. 
-P,' - 


Let us now relax the assumption that d^ lies in the observable 
space of (A'C)* The observability matrix is given by M' defined 
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by ( 2 . 3 - 23 ) , so if is unobservable, M' = 0_. Therefore, 
d^ lies in the nullspace of M' , or equivalently, in the detec- 
tion space of As Beard shows, every vector contained in 

the detection space of £ has the same detection order and de- 
tection space as so D is a detector gain for d^ as well as 

for £ (i.e., d^ is detectable in the traditional sense). There- 
fore, the second step of making D' a detector gain for d^ is 
unnecessary, and Ce* (t) will lie in a fixed direction. 

“S 

Finally, it remains to show that all the eigenvalues of 
A-DC are assignable. When making D a detector gain for V 
eigenvalues can be assigned using (2.3-17) since (A,C) is an 
observable pair. If d^ is unobservable with respect to (A',C'), 
then D' is unconstrained, and the remaining n- V eigenvalues are 
assignable by the free choice of D' by the auxiliary detector 
gain lemma. Now if d^ is observable for (A',C'), a total of 
only n-rank M' eigenvalues can be assigned while making D' a 
detector gain for d^ and using D” , the auxiliary detector gain 
for D' (D' is the auxiliary detector gain for D) . The remaining 
V eigenvalues are associated with the unobservable space of 
(A',C') which coincides with the detection space of ^ and have 
already been specified by making D a detector gain for £. Thus 
all the eigenvalues are assignable. 

One last remark helps to simplify the computation of D'. 

D' can be made a detector gain for A£ instead of d^ because 
Af = (mod E) , where E denotes the unobservable space of 
(A',C'). This means that the difference A£ - d^ lies in E. 

With this substitution, the need to compute (2.3-26) is avoided. 
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2.3.4 Sets of Events 


Up to this point, we have only cor.csrned ourselves with 
detectability of a single event. Now consider the set of event 
vectors | f ^ , . . . , f^"^. The set is considered "mutually detectable" 
if there exists a D that satisfies the conditions of detectabil- 
ity for all r events. If the set of vectors satisfy 

rank CF = r 

F = I ..... (2.3-32) 

with jjL defined for each f^ as in (2.3-14), then the events are 
called "output separable." Two events are not output separable 
if the output error for both events lie in the same direction. 

All sets of events considered in this v;ork were output separable. 

A necessary and sufficient condition for a set of output 
separable events to be mutually detectable is given by the follow- 
ing theorem; 

Group Detection Theorem — 

The output separable vectors \ / • • • ' are mutually 

detectable if and only if the sum of the individual detec- 
tion orders of the is equal to the croup detection order 
given by the dimension of the nullspaca of (n-rank M^) , 

where M' is defined as M' in (2.3-24) with f replaced by F. 

G 

The basic outline of the proof proceeds as follows. Since 

the set of events [f^, /f^.] are output separable, the events 

have nonintersecting subspaces (shown by Heard), and thus 
eigenvalues can be specified for each ..'.aile satisfying 
(2.3-20). An additional = rank eieer.values can be specified 
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by free choice of D' . Therefore, condition (ii) of detectability 
is satisfied if and only if Z V. = n - q' . If the set is not 
mutually detectable, n - q' - eigenvalues cannot be con- 

vj 1 

trolled. 

2. 3. 4.1 Actuator Set 

The design of the detection filter for a set of actuator 
events proceeds as follows: 

(1) determine the maximal generator for each b^ using (2.3-10) 

(2) form F as defined in (2.3-32) and divide the set up into out- 
put separable subsets 

(3) subdivide the set further, if necessary, until the subsets 
are also mutually detectable 

(4) solve (2.3-20) for each b^ to make D a detector gain for F 
while using the p!^s to specify XVi eigenvalues 

(5) specify the remaining eigenvalues of A-DC using D' 

Step 3 is unnecessary if the uncontrollable eigenvalues re- 
sulting from nonmutual detectability are satisfactory. Beard 
generates an algorithm that aids in finding these fixed eigen- 
values and can also expedite step 3. 

2. 3.4.2 Sensor Set 

One design method for a filter detecting sensor failures 
proceeds as follows: 

A 

(1) for each c., determine f. such that Cf. = e where the 

-i' -1 -1 -p^ 

sensors are linearly independent 

(2) form F as in (2.3-32) and subdivide into mutually detectable 
sets (events are automatically output separable by step 1) 
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form the set of vectors [ Af ^ , . . . , Af where ^ • w f 

is the set resulting from step 2; let A' and C be defined 

as in (2.3-22) with ^ replaced by F, ; categorize the Af. 

^1 

as follows, and remove those events falling under (iii) : 

(i) M'p ^ii ^ o 

(iii) 0 ; 

Fj^^= , . . . be the set resulting from step 3; let 

A' and C be defined as in (2.3-22) with f replaced by 

F, ; repeat step 2 with A, C, and f. replaced by A', C, 

■^2 

and the Af^ in category (i) ; remove Af^'s until the A^^ 
are mutually detectable with respect to (A',C*); if any 
Af^j^ are removed, some A^^ may move from category (ii) to 
(i) and mutual detectability with the new members must be 
rechecked 

solve (2.3-20) for each f. in the set F, resulting from 
step 4 to make D a detector gain for F, while specifying 

JC3 

2 eigenvalues 

solve (2.3-20) for each Af^ in the set F, to make D' a 

detector gain for AF with respect to (A',C) while 

specifying eigenvalues (Vj^ denotes the maxi- 

mod "^mod 

:“;um order a detection generator for A^^ can have) 
use D" to specify the remaining eigenvalues of 
vA'-D'C) = A-DC [D” is to D* as D' is to D in section 


2. 3. 4.1.] 
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A few remarks concerning step 3 are in order. If category 
(i) is the case, A£^ does not lie in the unobservable space 
of C with respect to A' and thus Af_^^ is detectable. For 
categories (ii) and (iii ) , A£^ does lie in the unobservable 
space so A^^ is not detectable; this is of no consequence, 
however, if D is a detector gain for A^^ as well as as is 
the case for category (ii) (see end of Section 2.3.3). The 
rank test given in step 3 determines whether the given sequence 
of vectors lie in a fixed direction — if they do, category (ii) 
applies . 

The assignment of eigenvalues in the last three steps 
is depicted in Fig. 2.2; ^ eigenvalues are assigned using 

D and n - ZV. are specified using D'; of che n -XV., 

pi 1 

T. V • are assigned using Dp and the rest using which 
^mod 

depends on D" (subscripts P and H denote p rticular and 
homogeneous solutions, respectively). 

2.4 Two-mode Design Example 

As a simple example of the analytic design procedure for 
the failure detection filter, consider the normal mode equations 

'VCW + wfVrW = Gn-aW . . . (2-4 -d 

where represents the amplitude of the r^^ mode, CO^. its 

natural frequency, and the amplitude of the input u. In a 

single degree of freedom mass-spring system, for instance, Y 

2 

IS just the displacement x and (JO is the spring constant k 
divided by the mass m. 
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In particular, consider the system of equations for 
two flexible modes of a beam given in the form (2.1-1); 



~i, _ o <PM) 0 <PtM yi 

'it o o (flM Y, 

i 


where force actuators at translation rate sensors 

ato(j^, 0(2 3-^® used [^^(Xj) represents value of the i^ mode 
shape at x^] . 

To simplify the design procedure, Beard shows that 
if the matrix A is replaced by 

A = A - T>C ( 2 . 4 - 3 ) 

with D arbitrary, the detector gain is unaffected. Following 
this suggestion, let's replace A by the simpler form 

O O O O 

03? O O O 

O O O O 

O O - 0)t o 
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Now let's proceed to design a detection filter that 
can detect failures of either actuator #1 or actuator #2. 

The two event vectors for this case are and b^, the two 
col’oitins of B. If the two actuators are not both at the same 
location or placed symmetrically about the center (for a free- 
free beam) , the directions Cb^ and Cb 2 will be linearly inde- 
pendent and thus output separable. To determine mutual de- 
tectability, the detection order for each £ must be found. 

For an appropriate choice of D, K is found from (2.3-24) 
to be equal to A, and C' using (2.3—22) becomes 


C' 


o 

o 


C. O C3 

Ci O Cm 


(2.4-5) 


where the C^'s are functions of the elements of Cb and C. Now 
form 



" r 

c' 


1 

0 

P 

0 

0 



0 C^. 0 Cm 

CV 


-(OTc* 0 0 

IN 


"Ol)? Ct, 0 Ch 0 

C'K' 


0000 

0 0 0 (L 

C’K' 


0 0 0 0 



1 

0 

0 

0 

0 

L_ 


(2.4-6) 


The two rows of C are linearly dependent so that the rank of 
M' is 2 and thus the detection order of f is 2 for both b^ and 
—2. Since is the null matrix, the group detection order is 
four and the events are mutually detectable. 
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The maximal generator for f is second order in this 
case and so must satisfy the set (2.3-10): 


Ca=o 

C^^'• 0 

T 

If £ is represented by [g^ gj^ g^ g^] , the first of (2.4-7) be- 
comes 

9t * O (2.4-3) 

<?2,9b -t O 

where ^ is used to denote ((3(j ) . Since we have assumed 
linear independence of the sensors, 


^iz 


O 


(2.4-9) 


and we must accept the trivial solution gj^=g^=0. The second 
equation in (2.4-7) transforms into 



-I 








- 

3 

u 



<P.._ 




(2.4-10) 


with and ^2 denoting the nonzero elements of b, and (2.4-9) 
guaranteeing existence of the inverted matrix. 

All that remains is to solve (2.3-20) with k = 2 for 
each event. If A is chosen for all the filter eigenvalues, 
p^ = X and P 2 = -2 A . The particular solution of D is found 
to be 
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(2.4-11) 






3,a 


(f-OCftSxc 




+ i?.x 


where ^ denotes |3j for b^. [Note that the original A is used 
in obtaining this result.] 

Since all four eigenvalues have been assigned using 
(2.3-20), D=D and the design is complete. If the gain were 
being designed for only one actuator, the remaining two eigen- 
values would have to be assigned by some other means. This 
is often not easy to do analytically using determinants, and 
is best done computationally using D', especially for higher 
order systems. 
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CHAPTER III 


COMPUTATIONAL DESIGN OF FILTER 

For systems of order four and greater, the analytic de- 
sign approach outlined in the previous chapter becomes very 
cumbersome. The design example was actually a relatively 
simple case facilitated by the fact that the same types of 
sensors were used and the filter was designed to detect two 
events. In most cases, it is not possible to solve for the 
detection generator explicitly nor is it possible to assign 

all filter eigenvalues using D alone (definitely not for 

P 

sensor failure events) . This chapter will present a design 
alternative based upon a process called "orthogonal reduction" 
which is best implemented on a computer. All of these algo- 
rithms were first proposed by Beard but also appear in Jones 
as well. 

3.1 Orthogonal Reduction 

Many of the processes involved in the filter design involve 
finding the rank of a matrix or finding a vector in the null- 
space of a certain matrix. The detection generator, for example, 
is in the nullspace of M' , the detection space of £, and the 
detection order is n-rankM' . The orthogonal reduction procedure 
is an iterative process which generates a positive semi-definite 
matrix, the range space of which coincides with the nullspace 
of a given matrix V, N(V). 

Let V be an mxn matrix given by 
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(3.1-1) 


where are the rows of V. Let A denote the symmetric matrix 

whose range space is to coincide with the nullspace of V. Each 
row of V is then tested to see whether it is orthogonal to the 
range space of V (the rows of V span N'*’(V) ) by defining the 
auxiliary vector 

Wi - n V; (3.1-2) 


If / 9.' Zi is not perpendicular to the columns of X2. ( -Q. is 
symmetric) and £i is reduced so that w = 0. This is repeated 
until the columns of iT. , R(fl.) , are all orthogonal to the rows 
of V, N^(V) . 

The reduction algorithm proceeds as follows ; 

(i) withJfl ^ from the previous iteration, form the auxiliary 
vector 


«i = n. V . 


(ii) if w^ = set 


■^i+1 -^i 

if w^ ^ set 

i+1 i 


w.w. 

— 1—1 

T 

W. V . 
—1 —1 


(3.1-3) 


The procedure must begin with iT. ^ positive definite to ensure 
a remains positive semi-definite and ^ 0. The rank of 
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the matrix is equal to the number of reductions performed. 


The algorithm for subroutine ORTRED is presented in Fig. 

3.1 with several modifications tailored to the types of matrices 
we will be dealing with. First of all, a row can be skipped 
in the reduction process if it is known that it is linearly 
dependent on the preceding rows. This is especially useful 
for matrices such as M' in (2.3-24) which are generated cycli- 
cally. When the auxiliary vector for a particular row of M', 

(Cj K'') for example, is zero, then this row is linearly de- 
pendent on the preceding rows of M' , Also, the remaining rows 
of M' that depend upon c ^ for i > Jl , are also de- 

pendent upon the preceding rows of M' . Thus, after the appear- 
ance of the first zero auxiliary vector, the rest of the rows 
that depend on Cj needn't be considered. [ICODEl(J) is set 
equal to one at this termination point.] When the termination 
points for all c^'s have been reached, M' is completely reduced. 

IC0DE3 and WL are used when computing the maximal generator. 
IC0DE3 is the row of C associated with the last nonzero auxiliary 
vector WL. The arrays IQ and WF are Used when calculating D'. 

IQ stores the orders of K associated with the termination points 
for each row of C, and WF stores the last nonzero auxiliary vec- 
tor for each row. 


3.2 Input failure event design 

The computation of D for actuators is accomplished by call- 
ing three subroutines: SEPDET, DETGEN, and DGAIN. Subroutine 

SEPDET first determines whether the events are output separable — 
it then determines whether the events are mutually detectable. 
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In the process of determining mutual detectability, SEPDET 
calls subroutine DETGEN to find the maximal generator for 
each event. Finally, subroutine DGAIN is called to calculate 
D using the detection orders and generators from SEPDET. In- 
puts to the main program FDFIL (see Appendix A) are the sys- 
tem matrix A and its order N, the measurement matrix C and 
number of sensors NS, the number of actuators NA, the matrix 
of event vectors FALL and number of events NF, and the filter 
eigenvalues EV. 

3.2.1 Subroutine SEPDET 

The flowchart for SEPDET appears in Fig. 3.2. The pro- 
gram first computes for each f^ according to (2.3-14) and 
replaces f^ by in FALL and CF. The rank of CF is found 

using ORTRED to see whether the events are separable. If 
not, the error message 'NOT SEPARABLE' is printed and the 
designer must remove the dependent events from FALL and start 
again . 

Once the events are output separable, the program pro- 
ceeds to determine whether the events are mutually detectable. 
The matrix M' in (2.3-24) is generated using C and K for the 
full set of events and the rank is found using ORTRED to 
determine the group detection order, = N - rankM^. The 
call to DETGEN yields the individual detection orders 
If V- 7 ^ 21 V • , the events are not mutually detectable, and 
one or more events will have to be removed and the process 
started over. 
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Fig. 3.2; Flowchart for SEPDET 











3.2.2 Subroutine DETGEN 


This routine will calculate the "maximal'* generator 
and detection order V for an event, or alternatively, the 
"maximum" generator and its order if (C,A) is an un- 

observable pair. The process is illustrated in Fig. 3.3 and 
begins with calculating M' in (2.3-24) for the single event 
and then using ORTRED to find the detection order , V =N-rankM ' . 
The resulting nullspace in this case is the detection space 
of £ and coincides with the range space of the final iT. . 

In addition to being in the detection space of the 
detection generator must satisfy (2.3-10). This is accom- 
plished by applying ORTRED to the matrix 

= [t 1 CKi • • i tK''"'] <3-1-4) 

Starting with the final XL ^ from the first step. Since C spans 
a subspace one dimension greater than C' and the process is 
Started with Xlf which is orthogonal to C ' , only one row of C 
will not be terminated when first encountered. If WL is the 
last nonzero auxiliary vector for this row , the maximal 
generator is given by 



As stated earlier, if (C,A) is not observable, the detec 
tion generator found in (3.1-5) is not maximal (order is not 
V), but it is the generator of maximum possible order; this 
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Fig. 3.3: Flowchart for DETGEN 











generator will be referred to as the "modified” generator and 
has order = rank If (C,A) is an observable pair, 

will equal V . The rest of the algorithm computes the 
vector CK^ which is one component of the output error for 
sensor failures. 

3.2.3 Subroutine DGAIN 


This subroutine will finally compute the detector gain 
D while assigning the filter eigenvalues. The process begins 
(see Fig. 3.4) by finding the particular solution of D for 

DCF = Q, 


where 





and 2.j^ have already been calculated by SEPDET. The 

‘^i 

coefficients p. are found by calling EVAS which computes the 
coefficients in (2.3-17) generated by (s-A) ^ where A is 

the chosen filter eigenvalue (all are same) . 

If = n, all eigenvalues are assigned and the de- 

tector gain is simply Dp. Otherwise, the remaining eigenvalues 
must be assigned using D' . Calculating D' is very similar to 
the procedure used in calculating Dp except that instead of 
using maximal generators, the last nonzero auxiliary vectors 
before termination found in reducing are used. Jones refers 
to this process as "completion of the state space” because one 
is assigning eigenvalues associated with that portion of the 
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state space not covered by a detection space. 

Finding D' begins by applying ORTRED to Mp (M' for the 
full set F) to obtain IQ and WF as explained in section 3.1. 
D' is then any solution satisfying 

where 



The Wp are the columns of WL corresponding to the rankC’ 

^i 

nonzero q., and Ai, C' 'correspond to A',C' in (2.3-22) for 

X r r 

the full set F. The coefficients p^j are found by calling 
EVAS for eigenvalues. One solution of (3.1-6) is simply 
the particular solution 

The general solution for D is found by applying (2.3-21) to 
the full set F. 

3.3 Measurement failure event design 

The algorithm for sensor filter design does not use any 
new subroutines (except SENSOR which calls the routines) but 
the order and implementation are different. A schematic for 
the sensor design process is given in Fig. 3.5. SEPDET is 
called first for (A,C,F) to determine whether the set is 
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restart 


- SEPDET - 
A, C, f, 


if IfAL-o 


IDET=1 


take out NMD f. 


- DGAIN - 

A. C> f. 

A’ 

\ f 

CATEGORIZE 
fi) . (ii) . (ii: 


IDET=1 


take out cat(iii)- 


- SEPDET - 
A' /Af 

ISEP=1 


take out NOS Af. 


IDET=1 



- SEPDET - 
A, C, 


- DGAIN - 
A, C, 

D =D 

pr 


D = D + D. 
P h 


take out NMD Af . 

X 


Fig. 3.5: Sensor design process schematic, 








mutually detectable (since ^ ~ lf...NF, it is auto- 

matically separable). If the set is not mutually detectable, 
events must be removed from F until it is and NF (no. events) 
set accordingly (F must now be inputted directly into the 
program instead of being calculated in SENSOR) . 

DGAIN is now called to calculate (set IAP=1) for use 
in computing Mp, the observable space of (C*,A*) (this has 
been denoted in the program) . The A£j|^ are then generated 
(recall that f is actually A^f when 0) and categorized 

according to step 3 in Section 2.3'. 4.2. The f^^ corresponding 
to the A^^ falling in category (iii) are then removed from 
F, and the process begins again with the new F and NF until 
no Af^ falls under (iii) . 

The next step is to call SEPDET for (A*,C*,A£^). If the 
A£^ are not separable, there is no problem since we still have 
the output error directions from the events to distinguish 
between sensor failures. Removing the dependent A^^ from the set 
AF^ we only need make D’ a detector gain for this new set AF 2 
of NF 2 events. [The program can be easily modified to handle 
this case by inputting AF 2 and NF 2 into SENSOR and bypassing the 
statements NF 2 =NF and AF 2 =AF . ] If the Afjj^ are not mutually de- 
tectable, however, the corresponding f^ will have to be removed 
from F and the whole process restarted with the new F. 

Once the A^^ are separable/detectable, DGAIN is called to 
calculate D' (D' is the output D of DGAIN in this case). SEPDET 
is recalled to regenerate and for the original to be 
used in calculating D^. Now DGAIN is called a final time for 
(A,C,fj|^) with the input D' calculated above to compute D. 



CHAPTER IV 


APPLICATION TO FLEXIBLE BEAM 

In demonstrating various control and failure detection 
algorithms for large space structures, researchers often use a 
long beam or rectangular plate as a typical structural element. 

A beam is usually chosen because of its simplicity and resem- 
blance to a long truss, though plates find their usefulness in 
simulating the closely spaced frequencies of a flat solar array. 

In this chapter, the experimental beam at NASA Langley Research 

/ 

Center will be described along with the design of the filter 
using the finite element description of the beam. The effec- 
tiveness of the filter in detecting failures of force actu- tors 
on the beam will then be tested using various input frequencies 
and filter bandwidths in mismatched filter-system models. A 
position sensor failure will also be simulated for the case where 
there is only an initial condition on modal amplitudes present. 
The effect of data sampling will also be investigated. 

4 . 1 NASA LaRC Experimental Beam 

In all the failure detection simulation tests, the dynam- 
ics model and actuator/sensor types and locations were chosen 
to correspond to those of the experimental beam at NASA LaRC 
to predict the performance of the filter in subsequent tests 
using the actual beam. The Langley beam is made of aluminum 
(M=0.502 slugs) and is twelve feet long with a 6" x 3/16" cross 
section (see reference 4) . Actuation is provided by electro- 
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magnetic shakers at four discrete locations along the beam, 
and the horizontal deflections are measured by nine Kaman 
KD-2300 probes, four of which are colocated with the actua- 
tors (see Fig. 4.1). All of the deflections of interest to 
us will take place in the horizontal plane of the beam. 

Because the beam is supported in a gravity field and 
the actuator dynamics are significant, the experimental beam 
does not behave exactly as a uniform beam would. Researchers 
at Langley have therefore found it convenient to perform their 
structural analysis using a finite element mode. The beam 
model consists of 24 sections, each of which is constrained 
so that it can translate and rotate in only one degree of 
freedom. The SPAR program was then used to obtain the natural 
frequencies and orthoconal mode shapes for the first ten modes 
of the beam (see Fig. 4.2a-h for first eight modes). The SPAR 
generated mode shape values at positions 0.5 ft, 2.5 ft, and 
6.0 ft from the end of the beam are given in Table 4.1 for the 
two rigid body modes (translation and rotation) and first six 
flex modes along with the modal frequencies. 

4 . 2 State equations for beam model 

In order to obtain the modal equations for the beam, we 
begin with the series expansion for the beam displacement 

^(e,0 - ^(P;(^)^(t) (4.2-1 

where is an orti.ogonal set of mode shapes and (t) are 
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Fig. 4.2; 


SPAR beam modal frequencies and shapes (frequencies in rad/sec) 
(taken from ref erence [ 4 ] ) 







the modal amplitudes. Next substitute this into the governing 
differential equation for the undamped beam 


^ (4,2-2) 

where f is the forcing term and m, E, and I are the beam mass 
(M) per unit length, modulus, and cross-section inertia, res- 
pectively. Assuming the use of m point force actuators, 

= Z! S(l-OU;W (4.2-3) 

J'-» 

with being the actuator positions and Uj (t) the control 
magnitudes, one obtains the set of relations 

cot ^ ~ WV ^ Clj = O (4.2-4) 

J 

where CO ^ is the natural frequency of the i^^ mode and is 
the normalized mode shape (for (j)^ not normalized, the input 
amplitudes must be scaled) . 

Next, consider casting (4.2-4) in the state-space form 

A 'X -i- "Bu (4.2-5) 

^ - C X 

» 

where 

[r y Y.t ■ X t] 
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where the niimber of modes has been truncated at N, and the 
use of m force actuators at positions and p translation 
deflection sensors at positions have been assumed. 


4.3 


Failure detection filter design for the beam 


For all the beam simulations that will appear later 


56 



in this chapter, the actuators and sensors located, at 0.5 ft, 

2.5 ft, and 6.0 ft from the end of the beam were used. The co- 
location was chosen simply as a matter of convenience, whereas 
the sensor located at 9.5 ft was not used since it is not out- 
put separable from the one located at 2.5 ft. In order for the 
state equations (4.2-5) to be a true representation of the beam, 
normalized mode shapes must be used. Thus, the orthogonal mode 
shapes in Table 4.1 were normalized according to (l = beam length) 

Ax = L (4.3-1) 

o 

before being used to compute A, B, and C. The orthogonal mode 
shape values at the actuator locations appear in Table 4.3. 

For actuator failure detection, it was always possible 
to design a filter that could detect a failure of two actuators, 
but never for all three. This is because the individual detec- 
tion orders are always two, but the group detection order for 
the three events is always n (the system order) making the 
events nonmutually detectable. The events are always output 
separable since the actuators are not symmetrically opposed with 
respect to the center of the beam. By extrapolation, it will 
be possible to detect up to p-1 input failure events with a 
single filter where p is the number of position sensors. 

The design for sensor failures was made somewhat more 
difficult by the fact that when designing for two sensors, it 
was found that Af^ and Af _2 were not output separable, with 
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respect to (C',A'). This simply meant that D' only needed 
to be made a detector gain for either or A ^2 would 

automatically be a detector gain for the other. The rest 
of the design proceeded normally resulting in a filter that 
could detect two failure events. As in the actuator case, 
it was not possible to detect three sensor events with a 
single filter. 

4 . 4 Computer simulation of beam 

In order to simulate the dynamics of the system and 
filter, a fifth order Runge-Kutta integration routine was used. 
The step size was always maintained at a sufficiently small 
size compared to the period of the highest frequency mode be- 
ing simulated in order not to misrepresent the continuity of 
the system dynamics. A listing of the beam simulation pro- 
gram appears in Appendix B. 

The program accepts as input the detector gain D cal- 
culated by FDFIL and the nominal system (filter) matrix A, 
actuator effectiveness matrix BF, and measurement matrix CF 
used in computing D. BS and CS for the system are also init- 
ially set equal to the nominal matrices BF and CF, respectiv- 
ely. When the failure time, TFAIL, is reached, the column 
of BS associated with the failed actuator, or correspondingly, 
the row of CS associated with the failed sensor, is zeroed 
out to represent a failure in the off position. The values 
of BF and CF retain the nominal system values thus creating 
the driving term in the state error equation. 
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To easily verify whether the steady-state output 
error lies in the direction Cb^ subsequent to the failure 
of the i^^ actuator, the output is transformed so that any 
vector lying in the direction Cb^ is changed to the direc- 
tion . The desired transformation then is 

-R’ [(CFyCFl'CCFy 

where F is the full matrix of event vectors b^^. It is 

easily verified by direct substitution of (4.4-1) that 

RCb^ = e for i=l,2,3. 

^i 

Similarly, it would be desirable to transform the 
output error for the i sensor event so that an error in 

A 0 

the plane spanned by and CK will lie in an easily 

identifiable plane. Fortunately, the use of three sensors 

permits a geometrical construction of the desired transforma- 

A 0 

tion. Consider the plane spanned by e^^ and w = CK in 

Fig. 4.3. This is the plane in which the output error is 
constrained to lie in the event of a sensor #1 failure. 

Now if the xyz coordinate system is rotated counterclockwise 
about the x-axis through the angle Y (the angle that the 
projection of w on the y-z plane makes with the y-axis) , the 
£-^7^ plane coincides with the plane spanned by e^^ and w. 

Thus under the transformation 



1 O O 

0 eos\ s'lnY 

o -sinV cosX 


(4.4-2) 
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an error in the plane spanned by e^ and w will have a first 
and second component but no third. 

In the next two sections, the effect of model error 
in the filter and data sampling of the output will be inves- 
tigated for the deterministic system. These two effects 
will be of the greatest concern in actual tests of the beam 
since the filter must necessarily incorporate a truncated 
model and data processor capability limits the sampling 
frequency. The effect of actuator, sensor, and plant noise 
will not be investigated here because they enter the realm 
of stochastic filter design which is beyond the scope of 
this work. Even after satisfying the requirements of being 
a detector gain, there is usually some freedom left in D 
that can be used to make the behavior of the detection filter 
approach that of the Kalman filter. Since no attempt was 
made here to optimize D for noisy systems under the constraints 
of being a detector gain, results of noise simulations would 
not reflect the true capabilities of the detection filter. 

4.4.1 Effect of model error 

Consider the case where the filter model contains 
fewer modes than the system model — this will always be the 
case for the actual beam. Let the state and output equations 
for the system be 
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and for the filter 

L- ^«X, +'B„a 

A A ^ 

'!ijL - Ckv 


(4.4-3) 


(4.4-4) 


where the subscripts M and U denote modeled and unmodeled 
modes, respectively. 

The new state error equations for the modeled modes 

are 

-T)CoXu (4 

and the coupled equations for the system and filter become 




a 


(4.4-6) 


The effect of and C^ on the system of (4.4-3) are known 
as control and observation spillover, respectively. Spillover 
was first pointed out by Balas (7) in his work on active con- 
trol of flexible structures. B^ in (4.4-3) has the effect of 
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driving the unmodeled modes, and creates an extra driving 
term in the state error equation (4.4-5) — ^however, neither 
type of spillover can cause filter instability since the 
triangularity of the dynamics matrix in (4.4-6) is preserved. 

The first two simulations performed were for a perfectly 
matched system and filter (no spillover present) in order to 
verify the filter design. The detection filter was set up to 
detect a failure of either actuator #1 at 0.5 ft from the end 
of the beam, or actuator #2 at 2.5 ft. Results for the simu- 
lation of failures of both actuators at 1.0 sec are shown in 
Fig. 4.5 for a five-mode model and 15 rad/sec filter band- 
width. All three actuators were driving the system with an 

-0.2t 

exponentially decaying sinusoid of the form, e ' sinCO^t, 
with an input frequency of CO^ = 20 rad/sec. The three trans- 
formed output errors are plotted on the same graph with the 
nonzero-signal in Fig. 4.5a representing the failure signal 
corresponding to the direction and, in Fig. 4.5b, the di- 
rection £ 2 * verification is good to within the number 

of significant figures in the detector gain used in the simu- 
lation. 

The effect of model error is dramatically illustrated 
in Fig. 4.6 for failure of the first actuator. Actuators 
#1,2,3 were operating at 20, 5, and 50 rad/sec, respectively, 
and the filter bandwidth was 15 rad/sec. Figure 4.6a is for 
a four-mode filter model and eight-mode system. The failure 
signal is not distinguishable from the other signals because 
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ERROR (FT) ERROR (FT) ERROR (FT) 

lOi .004 .012 -^0013 ^0004 .0005 .0014 ^0013 -v 0004 ' .0005 



TIME (SEC) 

(b) transformed output error for actuator #2 failure. 



?^00 TIs !s0 TtS iToO iTS iTsO UTS 2 . 


TIME (SEC) 

(c) real output error for actuator #2 failure. 
Figs.4.5a-c: Actuator failure for matched models. 
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(c) 6-mode filter 

V Figs.4.6a-c: Model error effect for different filter models. 
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the unmodeled fifth mode with a natural frequency of 61 
rad/sec is being excited heavily by the input frequency of 
50 rad/sec. This excitation is much less for the five- 
mode filter of Fig. 4.6b since the mode at 61 rad/sec is 
now included in the model, and the effect in Fig. 4.6c for 
the six-mode filter is barely noticeable. Therefore, one 
should include in his model at least the first mode whose 
natural frequency lies beyond the highest input frequency. 

The results of an eight-mode system/five-mode filter 
set of tests for various input frequencies and bandwidths are 
shown in Table 4.2. The peak failure signal for actuator #1 
is given along with the peak signal-to-noise ratio — ^the peak 

A A 

noise is the peak signal achieved in direction ^2 ^ 3 * 

It appears that for a constant filter bandwidth, the perform- 
ance of the filter drops off sharply as the input frequency 
is pushed farther past the filter poles. This is due to 
faster attenuation of the failure signal than the spillover 
noise outside the filter bandwidth. For a fixed input fre- 
quency, the signal-to-noise ratio reaches a maximum of 45 
for a bandwidth of 15 rad/sec and input frequency of 5 rad/sec 
and a maxim\am of 36 for a bandwidth of 25 rad/sec and 

= 20 rad/sec. Results for three filter bandwidths and 
three input frequencies are plotted in Pigs. 4.7 to 4.9. 

It is interesting to note than when the input frequency 
is being held constant, the trends we are observing are for 


66 



TABLE 4.2:' Peak failure signals and signal-to- 
noise ratios for various input 
frequencies and filter bandwidths 
in actuator failure tests. 


FILTER 

POLES 

(rad/sec) 

INPUT ] 
5 

FREQUENCY 

20 

(ra 

d/sec) 

50 

*PS 

+S/N 

PS 

S/N 

PS 

S/N 

5 

.0158 

16 

.00181 

4 

- 

- 

10 

.00595 

40 

.00150 

21 

.000284 

5 

15 

.00297^ 

45 

.00119 

31 

.000255 

5 

20 

.00175 

41 

. 000925 

35 

. 000238 

5 

25 

. 00114 

35 

. 000709 

36 

— 

— 

30 

.000804 

29 

. 000562 

23 

— 

““ 


*peak failure signal in 1.5— 2.0 sec interval 
■•■peak failure signal divided by peak noise 
in other two output directions (1.5— 2.0 sec) 
^the corresponding results with the detector 
gain designed for detecting only actuator#! 
failures were: PS=0. 00121,S/N=1. 5 


TABLE 4.3: Normalized mode shape values at 

sensor and actuator locations 
(taken from reference [5] ) . 


MODE 

SHAPE 

x=0. 

MODE SHAPE VALUE 
5 ft x=2.5 ft 

at 
x=6 , 

.0 ft 

<P. 

1, 

,00 

1.00 

1. 

00 


-1. 

59 

-1.01 

0. 

, 00 


-1. 

.60 

-0.123 

1. 

.21 


1. 

,32 

-0.876 

0. 

.00 

<Pf 

1. 

,04 

-1.30 

1. 

,41 


-0. 

.753 

1.09 

0. 

. 00 

<Pl 

-0, 

.465 

0.356 

1, 

.40 


-0, 

.181 

-0.553 

0. 

. 00 
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ERROR tPT) ERROR CFT) ERROR (FT) 

003 .0000 .0003 -sOOlO -^0003 .0004 .0011 ->.0020 -w0009 .0002 



(c) oO^=50 rad/sec 


Figs.4.7a-c: Actuator tests for X =-10 rad/sec. 
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ERROR (FT) ERROR (FT) ERROR (FT) 

-s00070 -^00035 .00000 .00035 ->0014 -n 0005 .0004 .0013 0030 -vOOiS .0004 




^ao 


IF 


Tio r?5 iTFo i75 uio iT^s H)Q 
TIME (SEC) 

(a) 0J^=5 rad/sec 




Figs.4,8a-c; Actuator tests for <X=-15 rad/sec. 
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ERROR (FT) ERROR (FT) ERROR (FT) 

-^0009 -^0095 ->.0001 .0003 ->,0020 -^0008 . OOOi .0018 ->.0060 -,.0025 .0010 .0045 





Figs.4.9a-c; Actuator tests for A =“20 rad/sec. 
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constant control spillover — ^however, the observation spill- 
over coefficient DC^ in (4.4-5) changes since D and the fil- 
ter poles change. But even if the filter poles are held 
constant and D is changed, the effects can be quite dramatic. 
For instance, if the detector gain is designed to monitor 
failures of just actuator #1 instead of both actuators, the 
signal-to-noise ratios for = 5 rad/sec between filter 
poles of 15 and 20 rad/sec are about two rather than 40 to 
50. The failure signal itself, however, is relatively un- 
affected. 

A set of simulations was also run for detection of 
complete failure (zero signal) of the sensor colocated with 
actuator #1. Instead of driving the system with actuator 
inputs, initial amplitudes were placed directly on the various 
modes in order to excite them. Since the output transforma- 
tion for the sensor case is a pure rotation of coordinate axes, 
the magnitude of the outputs will actually be in feet if 
normalized mode shapes are used. Thus, the transformed 
output errors have physical significance, whereas their mean- 
ing was obscured by the transformation in the actuator case. 

The results of the sensor failure simulations for filter 
poles at 15 rad/sec and varying initial conditions are given 
in Table 4.4 and Fig. 10. In all cases, the two rigid body 
modes had no initial conditions, and the amplitudes on the 
first three flex modes were 0.12 ft, 0.09 ft, and 0.06 ft, 
respectively. The corresponding filter initial conditions 
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TABLE 4.4; Sensor failure results for filter pol®s at 
15 rad/sec and various initial conditions. 


INITIAL 

AMPLITUDES 

(ft) 

X 

PEAK OUTPUT SIGNALS (ft) 

^S/N 


*PS: 

“PN 

.0045 

. 0030 

.0015 

.317 

.561 

5.92x10 

95 

. 0030 

.0020 

.0010 

.317 

.561 

3.95x10_3 

142 

. 0015 

. 0010 

. 0005 

.317 

.561 

1.97x10 . 

28 5 

.0003 

. 0002 

.0001 

.317 

.561 

3.13x10 ^ 

1790c 

0 

0 

0 

.317 

. 561 

5.36x10 . 

10^ 


TABLE 4.5: Sensor failure results for various filter 

bandwidths and fixed initial conditions 

(r6o=-003ft, * 


FILTER 

PEAK OUTPUT SIGNALS (ft) 


POLES 

4> 



V . 

(rad/sec) 

PSi 

*PS2 

PN 

S/N 

5 

.781 

1.191 

.00810 

147 

10 

.491 

.800 

. 00456 

175 

15 

.316 

.561 

. 00395 

142 

20 

.224 

. 396 

.00394 

100 

25 

.177 

. 286 

.00288 

99 

30 

.144 

.211 

.00287 

74 


peak signal in direction #1 (1. 5-2.0 sec) 
*same for direction #2 
same for direction #3 

^peak signal in direction #1 and #2 divided 
by peak noise in direction #3 
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Figs . 4 . lOa-c : Sensor tests for various initial conditions. 
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(c) 20 rad/sec 

Figs . 4 . lla-c ; Sensor tests for various filter bandwidths. 
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were all zero. The amplitudes on the three unmodeled modes 
(next three flex modes) are given in the table. The failure 
signal plane is x £ 2 . The model error noise is very sen- 
sitive to initial conditions on the unmodeled modes, vary- 
ing almost proportionally. The failure signal itself though 
is unaffected. 

Table 4.5 shows results when the amplitudes on the 
three unmodeled modes are held constant at 0.003, 0.002, and 
0.001 ft and the filter bandwidth is varied. The signal-to- 
noise ratio reaches a peak near 10 rad/sec where the attenua- 
tion of the spillover noise relative to the failure signal 
is greatest. Between 25 and 30 rad/sec, there was almost no 
attenuation of noise whereas there was a 25% reduction in 
the failure signal resulting in loss of signal-to-noise 
ratio. Plots for several cases appear in Fig. 4.11. 

4.4.2 Sampled-data systems 

Instead of monitoring a system continuously, it is 
usually the case because of computational considerations 
that sensor data is taken at regular intervals and the out- 
put held constant between sampling periods. This scheme is 
depicted in Fig. 4.4 for the detection filter. The switch 
at #1 closes once every sampling period to receive an up- 
date, and the signal stays constant between sampling times. 
Switch #2 closes at the same instant switch #1 closes, and 
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is shown in order to emphasize that the output error to 
the filter remains constant between sampling updates. 

The effect of data sampling on filter performance 
in detection of actuator #1 failures for filter poles at 
15 rad/sec and input frequencies of 20 rad/sec is shown 
in Table 4.6. As expected, the signal-to-noise ratio in- 
creases with sampling rate, but even at 160 hz, it is 
still a long way from the continuous sampling case where 
it is 31 (see Table 4.2). As illustrated in Fig. 4.12, 
the failure signal remains constant for different sampling 
rates, whereas the noise decreases for faster rates. 

As a final illustration of data sampling, consider 
the case where the sampling rate is held constant at 64 hz, 
and filter bandwidth is varied. Table 4.7 shows results 
for filter poles between 5 and 20 rad/sec. At 5 rad/sec, 
the error signal due to spillover noise reached values 
nearly equal to the failure signal. Between 10 and 20 
rad/sec the signal-to-noise ratio increased as it did 
for the continuous case. Finally, at 25 rad/sec the filter 
became unstable because the filter time constants were too 
fast for the sampled-error feedback to have a stabilizing 
effect. Plots for several cases appear in Fig. 4.13. 

It is important to understand that these data-sampling 
results are for a case in which the detection filter was de- 
signed on a continuous data basis , and then data sampling 
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TABLE 4.6: Data- sampling results for input frequencies 

of 20 rad/sec, filter poles at 15 rad/sec, 
and various sampling rates. 


SAMPLING 
FREQUENCY (hz) 

PEAK OUTPUT ERROR SIGNALS (*10^) 
+PS *PN2 “PN3 

^S/N 

50 

11.7 

3.41 

2.21 

3.4 

64 

11.7 

2.17 

0.70 

5.4 

96 

11.7 

1.53 . 

0.43 

7.6 

160 

11.7 

1.02 

0.25 

11.5 


TABLE 4.7: Data-sampling results for input frequencies 

of 20 rad/sec, 64 hz sampling rate, and 
various filter bandwidths. 


FILTER POLES 
(rad/sec) 

PEAK OUTPUT 
PS 

ERROR 

PN 2 

SIGNALS {-lo"*) 
PN3 

S/N 

5 

17.9 

7.49 

12.5 

1.4 

10 

14.0 

3.53 

2.0 

4.0 

15 

11.7 

1.53 

0.4 

5.4 

20 

9.3 

1.35 

0.3 

6.8 


+peak signal in direction #1 ( 1 . 5 - 2.0 sec) 
*same for direction #2 
“same for direction #3 

^peak signal in direction #1 divided by peak 
signal in direction #2 and #3 (1. 5-2.0 sec) 
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ERROR (FT) ERROR (FT) ERROR IFTl 

•^0014 -^.OOOS .0004 .0013 ^0014 ^0005 .0004 .0013 -%0013 -v0004 .0005 .0014 





Figs.4.12a-c: Data sampling tests for various sampling 
rates(filter bandwidth is 15 rad/sec). 
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(a) 5 rad/sec 


QD 




Figs.4.13a-c; Data sampling tests for various filter 
bandwidths (sampling rate = 64 hz) * 
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was imposed on the resulting filter. It is possible to 
design the failure detection filter on the basis of a 
sampled-data description of the dynamics. This protects 
against the instability observed in these results, but 
data sampling fundamentally limits the ability of the 
filter to restrict a continuous failure signature to a 
single direction or plane in the output space. 
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CHAPTER V 


SUMMARY AND CONCLUSIONS 

The failure detection filter perforins very well in the 
context of flexible structure dynamics when the dynamic model 
in the filter matches that of the system it is monitoring. 

It was found that the filter usually works quite well even 
when the description of the system includes modes which are 
not modeled in the filter. This performance rapidly diminishes, 
however, as the input frequencies become close to the natural 
frequencies of the unmodeled modes. For a given input fre- 
quency, an optimum filter bandwidth for maximum signal-to- 
noise ratio usually exists not far beyond the input frequency. 

It is not always the case, however, that the filter will 
perform satisfactorily in the presence of modeling error. 
Sometimes two different detector gains capable of detecting 
the same failure event and designed with the same filter poles 
will yield grossly different detection performances — in one 
case the signal-to-noise ratio might be 50 and in the other 
case it might be less by an order of magnitude. This effect 
is caused by the direct appearance of the detector gain in 
the filter input due to observation spillover, the same 
phenomenon encountered when designing control systems for 
flexible structures. 

The filter worked well for sensor failures, even when 
the initial amplitudes on the unmodeled modes were a significant 
fraction of those on the modeled modes (up to 7.5% were tested). 
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By varying filter bandwidth, it was possible to obtain very 
high signal-to-noise ratios for fixed initial conditions. 

As expected, satisfactory performance of the filter de- 
signed for continuous data processing when used in a sampled- 
data mode depended strongly on a very high sampling rate. 

Only when the sampling rate reached five times the natural 
frequency of the last system mode did the signal— to— noise 
ratio surpass 10. As can be anticipated, when the sampling 
rate to filter bandwidth ratio is too low, the phase lag 
caused by sampling tended to destabilize the system. 

These observations bring to mind several suggestions 
with regard to both the actual beam experiment and areas for 
further study: 

(1) To ensure adequate filter performance, it may be neces- 
sary to eliminate observation spillover as much as pos- 
sible — it might be possible to accomplish this using a 
phase lock loop^^ to comb out only the modeled modes 
from the output or by reducing the bandwidth of sensor 
data (one must be careful in doing this, however, not 
to attenuate the failure signal as well) . 

(2) It may be necessary to use a sampling rate of at least 
64 hz in order to get stable and adequate performance 
from the filter — this will depend upon the degree of ex- 
citation of the higher modes. 

(3) Even after specifying D to be a detector gain and assign- 
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ing all the filter eigenvalues, there is still a lot of 
freedom left in D (especially when D is designed to 
detect only one event) that could be used to reduce 
spillover noise — ^using this freedom would require a 
more analytic approach to solving for D and involve 
solving a set of nonlinear simultaneous equations when 
assigning eigenvalues. 

(4) The whole issue of stochastic failure detection filter 
theory has not even been addressed here — ^this will be 
an additional consideration in designing the filter when 
optimizing signal-to-noise ratios in the presence of 
plant, actuator, and sensor noise. 
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APPENDIX A 


FILE: FDFIL FORTRAN A 


VM/SP CONVERSATIONAL MONITOR SYSTEM 


C** 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c* 

c 

c 

c* 

c 


PROGRAM FDFIL COMPUTES THE DETECTOR GAIN 
OF THE FAILURE DETECTION FILTER 
FOR BOTH ACTUATOR AND SENSOR FAILURE EVENTS 


INPUT: N - NUMBER OF SYSTEM STATES 

NS - NUMBER OF SENSORS 
NF - NUMBER OF FAILURE EVENTS 

IFAL - (*1) TO INPUT EVENT VECTORS FOR SENSOR CASE 
(-0) THEN CF=I IS GENERATED FOR ALL SENSORS 
IAS - (»1) FOR ACTUATOR EVENTS 
(«2) FOR SENSOR EVENTS 
EV - FILTER EIGENVALUES 
A - SYSTEM MATRIX 
C - MEASUREMENT MATRIX 
FALL - MATRIX OF EVENT VECTORS 
(COLUMNS ARE EVENTS) 

OUTPUT; 0 -DETECTOR GAIN 

OPR - AUXILIARY DETECTOR GAIN 

MU - ARRAY OF SEPARABILITY ORDERS 

NU - ARRAY OF DETECTION ORDERS 

NUMOD - ARRAY OF MODIFIED DETECTION ORDERS 

10 - ARRAY OF COMPLETED STATE SPACE ORDERS 

CF - MATRIX OF OUTPUT ERROR DIRECTIONS FOR 

ACTUATOR FAILURES (COLUMNS ARE VECTORS) 
CZKLAF - THIS VECTOR AND THE UNIT VECTOR IN THE 
ITH DIRECTION SPAN THE PLANE OF OUTPUT 
ERRORS FOR SENSOR FAILURES 

IMSL SUBROUTINES THAT ARE USED: 

VMULFM(A,B.C) - C«A( TRANSPOSE ) *B 
VMULFP( A,B,C) - C*A*B(TRANSPOSE 
VMULFFU.B.C) - C*A*B 
LINV1F(A,A1) - AI»A( INVERSE) 

EIGRF(A.EV) - EV«EIGENVALUES(A) 


THIS IS THE MAIN PROGRAM FOR ACTUATOR EVENT FILTER 


INTEGER NU( 10).ICAT( 10) 

REAL D0UT(20.10) 

DOUBLE PRECISION ^ ^ ^ 

& A(20,20),C( 1O.2O),FALL(20, 10).CF( 10, 10),G(20. 10) 

& ZKF(20.20).CPF(10,20).ZMF(200,20).D(20. 10), 

& 0C(20.20).A0C(20,20),WK(440).REVAL(40),REVEC(800 

& DPR(20, 10) ,AP(20.20).EV,EPS.T0L 

C0MPLEX*16 EVAL(20).EVEC(20.20) 

EOUI VALENCE ( EVAL( 1 ) , REVAL( 1 ) ) , (EVEC( 1,1) ,REVEC( 1 ) ) 

DATA IN, 10/7,8/ 


FDF00010 
FDF00020 
FDF00030 
FDF00040 
FDF00050 
FDF00060 
FDF00070 
FDF00080 
FDF00090 
FDF00100 
FDF001 10 
FDF00120 
FDF00130 
FDF00140 
FDF00150 
FDF00160 
FDF00170 
FDF00180 
FDF00190 
FDF00200 
FDF00210 
FDF00220 
FDF00230 
FDF00240 
FDF00250 
FDF00260 
FDF00270 
FDF00280 
FDF00290 
FDF00300 
FDF00310 
FDF00320 
FDF00330 
FDF00340 
FDF00350 
FDF00360 
FDF00370 
FDF00380 
F0FOO390 
FDF00400 
•F0F00410 
FDF00420 
FDF00430 
FDF00440 
FDF00450 
FDF00460 
FDF00470 
FDF00480 
FDF00490 
FDF00500 
FDF00510 
FDF0O520 
F0FOO53O 
FDF00540 
FDF00550 
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READ(IN.IO) N,NS.NF,IFAL, lAS.EV 
READ(IN,20) ( (A( I , J) . J«1 , 16) . 1*1 . 16). 

& ((C(I. J). J«1, 16). 1*1.8). 

& ((FALL(I,J).J-1.8).I*1. 16) 

10 F0RMAT(5I5/E14.4) 

20 F0RMAT(32(8F1O.3/). 16(8F10.3/). 15(8F10.3/).8F10.3) 

WRITE(I0,25) N.NS.NF.IFAL.IAS.EV 

25 FORMAT ( 'N*' . I2.5X. 'NS*' . 12. 5X. 'NF*' . I2.5X. ' I FAL- ' . 12 . 5X, 

& 'IAS*' .I2/'EV«' .E14.4/) 

WRITE (10.26) ((A(I, J),J*1.8),I»1.8), 

& ((C(l.J).J*1.8).I*1.a). 

& ( (FALL( I . J) . J*1 .8) . 1*1 ,8) 

26 F0RMAT(//'A*'/8(8F1O.4/)//'C*'/8(8F1O.4/)//'FALL*V8(aF1O.4/)) 
IF(IAS.E0.2) GO TO 27 . 

CALL SEPDET(A.C.FALL.N.NF.NS, ISEP. IDET. IDIR.C.CF.G.NU.NUMOD. 

& 2KF.CPF.2MF) 

IF((ISEP.EO. D.OR.dDET.EO. D) GO TO 100 

CALL DGAIN(C.A,EV.G.CF,CPF .NU.NUMOO.N.NF.NS.O.O.DPR. 

& 2KF.2MF.AP.D) 

GO TO 28 

* CALL THE MAIN PROGRAM FOR SENSOR EVENT FILTER 

27 CALL SENSOR(C, A. EV.G.NU.CF .CPF .N.NF .NS, 

& 2KF.2MF. FALL. IF AL. IDET. ICAT.D) 

00 271 1*1,10 

IF(ICAT(I).E0.3) STOP 
271 CONTINUE 

IF(IDET.EO.O) GO TO 28 

WR I TE ( 10 , 277 ) IDET . ( ( FALL ( I , J ) , J* 1 . 4 ) . I = 1 . 6 ) 

277 FORMAT ( ' IDET*' . 1 1 // ' FALL* ' /6 ( 4E 14 . 4/ )// ) 

GO TO 200 

28 CALL VMULFF(D,C.N,NS,N,20, 10.DC.20. lER) 

CALL MATSUB(A,DC.20.20,N,N. ADC) 

WRITE(I0,29) ((ADC(I, J). J»1,4).I*1,4) 

29 FORMAT ( ' A-0C»'/4(4E 14 . 4/ )/ ) 

C* VERIFY THE PLACEMENT OF FILTER EIGENVALUES 

CALL EIGRF(ADC.N.20.0,REVAL.REVEC.20,WK. lER) 

WRITE(I0.3O) (EVAL(I).I*1, 16) 

30 FORMAT( 'CLOSED-LOOP EVS ARE: ' / 16( 2F 10 . 4/ )/ ) 

WRITE(I0,40) ((D(I.J).d»1.3).I»1, 16) 

40 FORMAT('GAIN MATRIX D* ' / 1 6 ( 3E 1 6 . 6/ ) // ) 

DO 48 1*1 ,N 
DO 48 J»1.NS 

DOUT(I.J)*SNGL(D(I.d)) 

48 CONTINUE 

WRITE( 11,50) ((00UT(I.J).d-1,3).I»1, 16) 

50 FORMATt 15(3E16.6/),3E16.6) 

GO TO 200 

100 WRITEdO. 1 10) ISEP. IDET 

110 FORMAT('THE SYSTEM IS NOT SEPARABLE OR DETECTABLE:'/ 

& ' ISEP* ' . 1 1 . 10X , 'IDET* '.ID 

200 STOP 
END 


FDF00560 
FOF00570 
FDF00580 
FDF00590 
FDF00600 
F0F00610 
FDF00620 
FDF00630 
FDF00640 
FDF00650 
FDF00660 
FDF00670 
FDF00680 
FOF00690 
FDF00700 
FDF00710 
FDF0O720 
FDF00730 
FDF00740 
FDF00750 
FDF00760 
FDF00770 
FDF00780 
FDF00790 
FDF00800 
FDF00810 
FDF00820 
FDF00830 
FDF00840 
FDF00850 
FDF00860 
FDF00870 
FDF00880 
FDF00890 
FDF00900 
FDF00910 
FDF00920 
FDF00930 
FDF00940 
FDF00950 
FDF00960 
FDF00970 
FDF00980 
FDF00990 
FDF01000 
FDF01010 
FDF01020 
FDF01030 
FDF01040 
FDF01050 
FDF01060 
FDF01070 
FDF01080 
FDF01090 
FDF01 100 
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C» SEPOET DETERMINES SEPARABILITY ANO DETECTABILITY OF EVENTS 
C 

SUBROUTINE SEPOET( A . C , FALL . N, NF , NS . ISEP, IDET , lOIR.CSYS.CF ,6,NU, 
& NUMOO.ZKF.CPF.2MF) 

INTEGER MU( 10) .NU( 10) . I0( 10) ,NUMOO( 10) 

DOUBLE PRECISION 

& A(20.20).C(10,20),FALL(20. 10),CF( 10, 10) ,CPF( 10,20) . 

& G(20, 10).ZKF(20,20),2MF(200,20),F(20).CAJ( 10,20), 

& CAUF( 10) ,0M( 10, 10) ,0MSU8( 10, 10) , WL( 20) , WF( 20. 10) , 

& CFTC( 10,20) ,CFTCF( 10. 10) . CFTCF I ( 10. 10) . WKAR( 460) . 

& CFTTI ( 10,20) ,AFALL( 20, 10) , AFCF( 20. 20) , CFCF( 10.20) . 

& 2MFNUL(20.20).AJ(20.20).ZMFSUB(20.20).FALL2(20. 10). 

& ZMFSB( 10.20),ZMFSB2( 10 . 20) . A J2( 20. 20) . AMUF( 20) , 

& CSYS( 10,20).EPS.T0L,CAUEPS.ZKNU1 

DATA EPS. 10/1 .0-4,8/ 

ISEP»0 
IDET=0 
DO 3 I=«1.N 
DO 3 U»1 ,NF 

FALL2(I .d)»FALL(I. J) 

3 CONTINUE 
DO 20 K=1 ,NF 
DO 5 1=1 ,N 

F(I)-FALL(I,K) 

5 CONTINUE 
MU(K)»0 

CALL MATID(AJ.20.N) 

6 CALL VMULFF(C.Ad,NS.N.N. 10.20.CAJ, 10, lER) 

CALL MATVEC(CAd.F. 10. 20. NS , N, CA JF ) 

TOL=O.DO 

DO 7 1=1, NS 
DO 7 J=1.N 

CAJEPS=CAd(I. J)*EPS 
TOL=OABS( DMAX 1 ( TOL , CA JEPS ) ) 

7 CONTINUE 
NZERO-O 

DO 10 1=1, NS 

IF(DABS(CAvlF(I)).LT.TOL) NZERO=NZERO+ 1 

10 CONTINUE 
IF(NZERO.NE.NS) GO TO 12 
MU(K)=MU(K)+1 

CALL VMULFF(A.AJ.N.N.N,20.20.AJ2.20. lER) 

DO 11 1=1. N 
DO 11 J-1.N 

AJ(I.J)»AU2(I.d) 

11 CONTINUE 
GO TO 6 

12 WRITE(I0.13) K.MU(K) 

13 F0RMAT('K-'.I2.5X.'MU«'.I2) 

CALL MATVEC(Ad,F,20.20.N,N.AMUF) 

DO 15 I-1.N 

FALLd ,K)«AMUF(I) 

CF(I.K)»CAJF(I) 

15 CONTINUE 
20 CONTINUE 


F0F01 1 10 
FDF01120 
FDF01130 
FDF01 140 
FDF01150 
FDF01 160 
FDF01170 
F0F01 180 
FDF01190 
FOF01200 
F0F01210 
FOF01220 
FDF01230 
FDF01240 
FDF01250 
FDF01260 
FDF01270 
FDF01280 
FDF01290 
FDF01300 
F0F01310 
FDF01320 
FDF01330 
FDF01340 
FDF01350 
FDF01360 
FDF01370 
FDF01380 
FDF01390 
FDF01400 
FDF01410 
FDF01420 
FDF01430 
FDF01440 
FDF01450 
FDF01460 
FDF01470 
FDF01480 
FDF01490 
FDF01500 
F0FO151O 
FOF01520 
FDF01530 
FDF01540 
FDF01550 
F0F01560 
FDF01570 
F0F01580 
FDF01590 
F0F01600 
F0F01610 
FDF01620 
FDF01630 
F0F01640 
F0F01650 


86 



FILE: FDFIL 


FORTRAN A 


VM/SP CONVERSATIONAL MONITOR SYSTEM 


WRITE ( 10,205) ((CF(I.d).J-1.3).I»1,3) 

)5 FORMAT( 'CF-'/3(3E16.6/)/) 

CALL MATIDIOM. 10. NS) 

CALL ORTRED(CF. 10, 10. NS , NF . OM, NS . 0 , OM , OMSUB , 

& 1RKCF.IC0DE3.ZKNU1.WL.I0.WF) 

WRITE(I0,2O7) IRKCF 
207 FORMAT! 'RANK OF IS '.12) 

IF( IRKCF. EO.NF) GO TO 21 
ISEP-1 

WRITE! 10,209) 

39 FORMAT! 'THE EVENTS ARE NOT SEPARABLE') 

RETURN 

21 CALL VMULFM!CF.C,NS.NF.N, 10, 10.CFTC, 10, lER) 

CALL VMULFM!CF.CF,NS.NF.NF. 10, 10.CFTCF , 10, lER) 

CALL LINV1F!CFTCF,NF. 10.CFTCFI .O.WKAR. lER) 

CALL VMULFF!CFTCFI.CFTC.NF.NF.N. 10. 10.CFTTI . 10, lER) 

CALL VMULFF!A,FALL,N,N,NF.20.20.AFALL,20. lER) 

CALL VMULFF!AFALL.CFTT1 .N,NF,N.20, 10 . AFCF , 20. lER ) 

CALL MATSUB!A.AFCF.20.20.N.N.ZKF) 

CALL VMULFF!CF.CFTTI .NS.NF.N. 10. 10.CFCF, 10. lER) 

CALL MATSUB!C,CFCF, 10, 20. NS . N, CPF ) 

CALL MATGEN!CPF , ZKF . 10. 20, 200, NS , N , N, ZMFSB , ZMFSB2. ZMF ) 

WR I T E ! 1 0 . 303 ) !!ZKF!l.d).J=1,4),I = 1,4). 

& !!CPF!I.d).J=1.4).I=1,3) 

33 FORMAT! ' ZKF= ’/4! 4E 14 . 4/ )// ' CPF« ' /3! 4E 14 . 4/ )// ) 

CALL MATID!ZMFNUL.20.N) 

CALL ORTRED! ZMF , 200 , 20 , NS *N , N , ZMFNUL , NS , O , ZMFNUL . ZMF SUB , 

& IHKZMF.IC0DE3.ZKNU1,WL,I0,WF) 

CALL DETGEN!MU,FALL,FALL2,IDIR,CSYS, A,C,N,NS.NF,ZKF.NU,NUMOD.G) 

WRITE! 10.315) IRKZMF.!NU!I).I»1,2) 

15 FORMAT! 'IRKZMF«' .12. 10X. 'NU 1 -2= ' . 215// ) 

NUG’N-IRKZMF 
INUSUM=0 
00 25 I=1,NF 

INUSUM»INUSUM+NU!I ) 

25 CONTINUE 

IF! INUSUM. EO.NUG) GO TO 30 
IDET=1 
30 RETURN 
END 


C* DETGEN FINDS THE DETECTION ORDER AND GENERATOR FOR THE EVENTS 
C 

SUBROUTINE DETGEN! MU. FALL , FALL2 . IDIR . CSYS , A , C , N . NS . NF . ZKF . 

& NU.NUMOO.G) 

INTEGER MU! 10) ,NU! 10) , lO! 10).NUMOO! 10).LEXP! 10) 

DOUBLE PRECISION 

& FALL! 20. 10). A! 20. 20) ,C! 10,20) .CSUB! 10.20). 

& G!20, 10).F!20).CSF! 10) . CSFTC ! 20) . CT ! 20 . 10). 

& ASF!20) .ASUB!20,20).ZK!20.20) .CP! 10. 20) . ZMP! 200. 20) . 

& ZMK!200,20) .ZMKNUL! 10. 10) . Cl ! 20) . AMU! 20. 20) . 

& WL!20).WF!20. 10).ZKF!20,20) .ZMPNUL!20.20).ZMKSUB!20.20) 

& ZMPSUB!20.20) .ZMPSB! 10. 20) , ZMKSB! 10.20) . 

& ZMPSB2! 10,20) .ZMKSB2! 10.20) ,AMU2!20, 20) . 

& AF!20).FALL2!20. 10) . ZKL! 20 . 20) . CZKL! 10,20). 
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& CZKLAF( 10) ,2KL2( 20,20) ,CSYS( 10 , 20) , C2KEPS , 

& CFTCF , 2KNU 1 , Cl AMUF . EPS . TOL 

DATA 10. EPS/8, 1.0-4/ 

DO 45 K»1 ,NF 
DO 5 1=1. N 

F(I)-FALL(I.K) 

5 CONTINUE 

CALL MATVEC(C,F. 10,20,NS.N,CSF) 

DO 3 1=1, NS 
DO 3 J»1,N 

CT(J.I)»C(I. J) 

3 CONTINUE 

CALL MATVEC(CT.CSF.20, 10. N . NS . CSFTC ) 

CFTCF»0. 

DO 10 1*1 ,NS 

CFTCF=CFTCF+CSF( I )**2 
10 CONTINUE 

DO 15 1 = 1. N 

CSFTC( I )«CSFTC(I )/CFTCF 
15 CONTINUE 

CALL MATVEC(A.F,20.20.N.N.A5F) 

DO 17 I»1,N 
DO 17 J=1,N 

ASUB( I , J)=ASF(I )*CSFTC(U) 

17 CONTINUE 

DO 18 1=1. NS 
DO 18 J=1,N 

CSUBd . J)=CSF(I)*CSFTC( J) 

18 CONTINUE 

CALL MATSUB(A,ASU8.20,20,N,N,ZK) 

CALL MATSUB(C.CSUB, 10 . 20 , NS . N . CP ) 

CALL MATGEN( CP , ZK . 10 . 20 , 200 . NS , N . N , ZMPSB , ZMPSB2 , ZMP ) 
CALL MATID(ZMPNUL.20,N) 

CALL ORTRED( ZMP , 200 . 20 . NS *N . N , ZMPNUL . NS , 0 , ZMPNUL . ZMPSUB . 
& I RKZMP . I C0DE3 . ZKNU 1 , WL , I 0 . WF ) 

WRITE( 10,206) I RKZMP 
206 F0RMAT('RANK OF ZMP IS dI2) 

NU(K)=N-IRKZMP 

CALL MATGEN(C.2K, 10. 20 , 200, NS . N. NU(K) . ZMKSB . ZMKSB2 , ZMK ) 
CALL ORTRED(ZMK,200,20,NS*N,N,ZMPNUL,NS.O,ZMKNUL,ZMKSUB, 
& IRKZMK, IC00E3.ZKNU1 .WL, lO.WF) 

NUM00(K)*IRK2MK 

WRITE( 10.216) NUMOO(K) 

216 FORMAT( 'NUM00»'.I2) 

DO 20 J*1,N 

CI( J)*C(IC00E3. J) 

20 CONTIMUE 

Cl AMUF =0.00 
DO 25 1=1. N 

CIAMUF-CIAMUF+CI ( I ) *F ( I ) 

25 CONTINUE 

DO 40 I=1,N 

G( I .K)»(CIAMUF/ZKNU1 )*WL( I ) 

40 CONTINUE 

43 WRITE(I0.225) ( G( I . K) , I = 1 . 4 ) 


FDF02210 

FDF02220 

FDF02230 

FDF02240 

FDF02250 

FDF02260 

FDF02270 

FDF02280 

FDF02290 

FDF02300 

FDF02310 

FDF02320 

FDF02330 

FDF02340 

FDF02350 

FDF02360 

FDF02370 

FDF02380 

FDF02390 

FDF02400 

FDF02410 

FDF02420 

FDF02430 

FDF02440 

FDF02450 

FDF02460 

FDF02470 

FDF02480 

FDF02490 

FDF02500 

FDF02510 

FDF02520 

FDF02530 

FDF02540 

FDF02550 

FDF02560 

FDF02570 

FDF02580 

FDF02590 

F0FO2600 

FDF02610 

F0FO262O 

FDF02630 

FOF02640 

FDF02650 

FDF02660 

FDF02670 

FDF02680 

F0F02690 

FDF02700 

FDF02710 

FDF02720 

FDF02730 

FDF02740 

F0FO275O 
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225 


441 


442 


443 


444 


445 

446 

447 


FORMAT( 'G»'/4E14.4) 

IFdDIR.NE. 1) GO TO 45 
LEXP(K)»0 

00 441 I«1.N 

AF(I)-FALL2(I.K) 

CONTINUE 

CALL MATI0(ZKL.2O.N) 

CALL VMULFF(CSYS.ZKL,NS,N.N, 10,20.CZKL, 10. lER) 

CALL MATVECICZKL.AF. 10.20.NS.N.CZKLAF) 

TOL-0.00 
DO 443 I«1.NS 
00 443 J»1.N 

CZKEPS“CZKL( I , J)»EPS 
TOL=OABS(OMAX1(TOL.CZKEPS)) 

CONTINUE 

NZER0=O 

DO 444 I»1.NS ^ 

I F ( OABS(CZKLAF( I ) ) . LT . TOL ) NZERO=NZERO+ 1 

CONTINUE 

IF(NZERO.NE.NS) GO TO 446 
LEXP(K)=LEXP(K)+1 

CALL VMULFF(ZK.ZKL.N.N,N,20,20,ZKL2,20. lER) 

DO 445 I»1.N 
DO 445 J“1.N 

ZKL(I.vJ)-ZKL2(l.J) 

CONTINUE 

GO TO 442 . 

WRITE (10.447) UEXP( K ) . ( CZKLAF ( I ) . I = 1 . NS ) 
FORMAT ( ' L» ' . 12 , 5X, ' CZKLAF= ' , 3E 16 . 6) 


C 

C* 

C 


45 CONTINUE 
RETURN 
END 

DGAIN DETERMINES THE DETECTOR GAIN FOR THE FILTER 

SUBROUTINE DGAIN( C . A . EV . GEN . CF . CPF , NU , NUMOO . N. NF . NS . I AP . lOPR . DPR 

& ZKF.ZMF.AP.D) 

INTEGER NU( 10) . I0( 10) .NUMOD( 10) 

DOUBLE PRECISION 


& C( 10,20) . A(20,20) ,GEN(20, 10) ,CF( 10. ^ * 

& G(20) .0(20. 10) .P(20) .0(20) ,Ad( 20.20) ,AOG( 20) . 

I 00(20 10).CFTCF( 10.10). CFTCFH , 

& CFTTI( 10. 10) , DP (20, 10) , ZMF ( 200, 20) • J ’ 

& WL(20) WF(20. 10) . AP(20. 20) . DPC( 20. 20) .PSI(20) . 

& APO(20.20).WFI(20).APJWF(20).PSIW(20. 10) . ZKF ( 20. 20) . 

A ZKF0(20 20) ZKFWF(20) .W(20, 10) .CPFW( 10. 10) .CPFWT( 10. 10) 

I ?JfSti?;o l6) CFWm(10.10).DPR(20.10).CF2TTI(10 

A DPRSUB(20. 10) .0H(2O. 10) .OMSUB( 10. 10) . AJ2(20.20) . 

I AP J2( 20. 20) . ZKF02( 20. 20) . ZMFSUBC 20. 20) . EV . EPS . TOL 

DATA 10/8/ 

DO 35 K»1 .NF 
DO 5 I * 1 . N 

G(1)-GEN(I.K) 

5 CONTINUE 

CALL EVAS(20.NUM0D(K).EV.P) 


F0F02760 

FDF02770 

FDF02780 

F0F02790 

F0F02800 

F0F02810 

FDF02820 

F0F02830 

F0F02840 

FDF02850 

FDF02860 

FDF02870 

FDF02880 

FDF02890 

F0F02900 

FDF02910 

FDF02920 

FDF02930 

FDF02940 

FDF02950 

FDF02960 

FDF02970 

FDF02980 

FDF02990 

FDF03000 

FDF03010 

FDF03020 

FDF03030 

FDF03040 

FDF03050 

FDF03060 

FDF03070 

FDF03080 

FDF03090 

F0F03100 

FDF031 10 

FDF03120 

FDF03130 

FDF03140 

FDF03150 

FDF03160 

FDF03170 

F0F0318O 

FDF03190 

F0F03200 

F0F03210 

FDF03220 

F0FO323O 

F0F03240 

F0F03250 

FDF03260 

FDF03270 

FDF03280 

F0F03290 

F0F033OO 
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DO 10 I* 1 ,N 

0 (I)»P( 1)*G(I) 

10 CONTINUE 

NUMOD 1 «NUMOD ( K ) - 1 
IF(NUM001 .EO.O) GO TO 21 
DO 20 0»1. NUMOD 1 

CALL MATP0W(A,20.N.d.AJ2.AJ) 
CALL MATVEC(Ad.G.20.20.N.N.AdG) 


DO 15 I«1.N 

0(I)*0(I)+P(d+1)*AdG(I) 

15 CONTINUE 

20 CONTINUE 

21 CALL MATPOW( A,20.N,NUM0D(K) . Ad2. AJ) 
CALL MATVEC(Ad.G.20.20.N.N.AJG) 

DO 25 I=1,N 

0(I)«0(I)+AJG(I) 


25 CONTINUE 

00 30 I«1.N 
0D(I.K)=0(I) 

30 CONTINUE 

35 CONTINUE _ 

CALL VMULFM(CF , CF , NS , NF ,NF , 10, 10,CFTCF , 10, lER) 

CALL LINV1F(CFTCF.NF. 10.CFTCFI .O.WK, lER) 

CALL VMULFP(CFTCFI,CF,NF.NF.NS, 10. 10.CFTTI . 10. lER) 
CALL VMULFF(OO.CFTTI .N.NF.NS.20. 10. DP . 20. I ER ) 

CALL VMULFF(DP.C.N.NS.N.20. 10.DPC . 20. I ER ) 

CALL MATSU8(A.0PC.20.20.N.N.AP) 

WRITE(I0.391) ((DP(I.J).d=1.3),I=1.A), 

& ((AP(I.J).J=1.4).I-1.4) 

394 FORMAT ( 'DP«'/4(3E14.4/)//'AP='/4(4E14.4/)/) 


NUSUM=0 
DO 40 K» 1 .NF 

NUSUM=NUSUM+NU(K) 


40 CONTINUE 

IFdDPR.EO. 1) GO TO 90 

IF( (N-NUSUM.NE.O) .AND. (lAP.NE. 

DO 45 I»1.N 
DO 45 d=1.NS 
0(1 ,d)=DP(r.d) 


)) GO TO 50 


45 CONTINUE 
RETURN 

CALL 0RTRED(ZMF*!'2OO?2o!nS*N.N.ZMFNUL,NS. 1 .2MFNUL.2MFSUB. 
8 . IRKZMF.IC00E3.2KNU1.WL.I0.WF) 

WRITE(I0.411) (10(1). 1*1. 4) 

411 FORMAT( 'I0«'/4I4) 


ICP»0 

DO 85 1-1 .NS 

IF(IOd).EO.O) GO TO 85 
lCP-lCP+1 

CALL EVAS(20.I0(I).EV.P) 
DO 55 L»1 ,N 

WFKL)-WF(L.I) 

55 CONTINUE 

DO 60 d-I.N 


FDF03310 

FDF03320 

FDF03330 

F0F03340 

FDF03350 

F0F03360 

FDF03370 

F0F03380 

FDF03390 

F0F03400 

FDF03410 

FDFO3420 

FDF03430 

FDF03440 

FDF03450 

FDF03460 

F0FO3470 

FDF03480 

F0F03490 

FDF03500 

FDF03510 

FDF03520 

FDF03530 

FDF03540 

FDF03550 

FDF03560 

FDF03570 

FDF03580 

FDF03590 

FDF03600 

FDF03610 

FDF03620 

FDF03630 

FDF03640 

FDF03650 

FDF03660 

FDF03670 

FDF03680 

FDF03690 

F0FO37OO 

FDF03710 

FDF03720 

F0F03730 

FDF03740 

F0F03750 

FDF03760 

FDF03770 

FDF03780 

FDF03790 

F0F03800 

FDF03810 

FDF03820 

F0F03830 

F0F03840 

FDF03850 
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PSK J)»P( 1 )*WFI(0) 

60 CONTINUE 

I0IM1»I0( I )-1 
IF(IOIMI.EO.O) GO TO 71 
00 70 U-1.I0IM1 

CALL MATP0W(AP,20.N. J,APU2.APJ) 

CALL MATVEC(AP0.WFI.2O.2O.N.N,APdWF) 

DO 65 K-1.N 

PSI (K)«PSI (K)+P(d+1 )»APJWF(K) 

65 CONTINUE 

70 CONTINUE 

71 CALL MATP0W(AP.20.N. 10(1 ) ,APd2. APvJ) 

CALL MATVEC(APd,WFI .20.20, N.N.APJWF ) 

00 75 K«1.N 

PSI(K)=-PSI(K)+APdWF(K) 

75 CONTINUE 

CALL MATPOWC ZKF . 20.N. I0( I ) - 1 . ZKFQ2 . ZKFO) 

CALL MATVEC(ZKF0.WFI.20.20.N.N.ZKFWF) 


DO 80 L«1.N 

PSIW(L. ICP)«PSI(L) 

W(L, ICP)-ZKFWF(L) 

80 CONTINUE 
85 CONTINUE 

CALL VMULFF(CPF.W.NS.N.ICP.10.20.CPFW,10.IER) 

CALL VMULFM(CPFW.CPFW.NS. ICP. ICP . 10, 10.CPFWT , 10, lER) 
CALL LINV1F(CPFWT.ICP. lO.CPFWTI.O.WK.IER) 

CALL VMULFP(CPFWTI .CPFW. ICP. ICP. NS, 10, lO.CFWTTI . 10, lER) 
CALL VMULFF ( PSIW . CFWTTI . N, ICP . NS . 20 . 10 . DPR . 20 . I ER ) 

90 WRITE( 10.457) ( (DPR( I . d) . J= 1 . 3 ) . I = 1 , 4 ) 

457 F0RMAT('DPR='/4(3E14.4/)/) 

CALL VMULFF(CF.CFTTI .NS.NF.NS. 10. 10.CF2TTI . 10. lER) 

CALL VMULFF(DPR.CF2TTI.N.NS.NS.20, 10 . DPRSUB . 20. lER) 

CALL MATSUB(DPR. DPRSUB. 20, 10. N, NS. OH) 

CALL MATA0D(DP,DH.20. lO.N.NS.D) 


RETURN 

END 


C 

C* 

C 


SENSOR IS THE MAIN PROGRAM FOR SENSOR EVENT FILTER DESIGN 

SUBROUTINE SENSOR(C.A.EV.G,NU.CF.CPF,N.NF .NS. 

4 ZKF. ZMF, FALL, IFAL. IDET, ICAT.D) 

INTEGER NU( 10) ,I0( 10) ,NUMOO( 10) ,ICAT( 10) 

DOUBLE PRECISION , ^ 

& A(20.20).C(10.20).G(20.10).CF(10.10).CPF(10.20). 

& ZKF(20.20) . 2MF(200.20) , FALL(20, 10) ,0(20, 10) ,CCT( 10. 10) 

& CCTI( 10. 10) , AP(20,20) .W0SUB( 10, 20) ,W0SUB2( 10. 20) . 

& WO(200.20).F(20).AF(20).WOAF(200).APAF(20).WK(40). 

& APAF2(20).CAPAF(10).CAAFN1( 10,20).CAFNUL(20.20). 

4 CAFSUB(20,20) ,WL(20) ,WF(20, 10) , AFALL(20. 10) . 

4 DPR(20. 10).CPF2( 10. 20) . AP2( 20. 20) . AFALL2(20. 10). 

4 EPS.TOL.WOEPS 

DATA ID. EPS/8. 1.0-4/ 

IOIR-0 

IF(IFAL.EO.I) GO TO 5 

CALL VMULFP(C.C.NS,N.NS. 10. 10.CCT , 10. lER) 


FDF03860 
FDF03870 
FDF03880 
FDF03890 
FDF03900 
FDF03910 
FDF03920 
FDF03930 
FDF03940 
F0F03950 
FDF03960 
FDF03970 
FDF03980 
FDF03990 
FDF04000 
FDF04010 
FDF04020 
FDF04030 
FDF04040 
FDF04050 
FDF04060 
FDF04070 
FDF04080 
FDF04090 
F0F04100 
FDF041 10 
FDF04120 
FDF04130 
FDF04140 
FDF04150 
FDF04160 
FDF04170 
FDF04180 
FDF04190 
FDF04200 
FDF04210 
FDF04220 
FDF04230 
FDF04240 
FDF04250 
FDF04260 
F0F04270 
FDF04280 
FDF04290 
. FDF04300 
FDF04310 
FDF04320 
FDF04330 
FDF04340 
F0F04350 
FDF04360 
F0F04370 
FDF04380 
F0FO4390 
FDF04400 
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121 


132 


CALL LINV1F(CCT.NS. 10.CCTI .O.WK. lER) 

CALL VMULFM(C.CCTI.NS.N.NS. 10, 10. FALL . 20 , 1 ER ) 

WRITE (10. 121) ((FALLd. J).d-1.A).I»1 .A) 

8, ZKF.CPF.ZMF) 

C ALL°OGA I N ( C ! a?ev'^g'!cF . CPF . NU . NUMOO , N . NF . NS . 1 , 0 . DPR . 

& ZKF.ZMF.AP.O) 

CALL MATGENC CPF , AP , 10, 20. 200, NS , N , N , WOSUB , W0SUB2 , WO ) 
WRITE (10. 132) ((W0(I.J).d=1.A).I»1. 12) 

FORMAT ( 'W0='/12(4E14,4/)//) 


NSN=NS»N 

T0L»O. 

DO 10 I»1.NS 
DO 10 d-1.N 

W0EPS-W0(I ,d)*EPS 
T0L*DABS(DMAX1(T0L,W0EPS) ) 


10 CONTINUE 

WRITE! 10. 143) TOL 
143 FORMAT! 'T0L«* ,E14. 4/) 

DO 60 K» 1 ,NF 
DO 15 I « 1 . N 

F(I)«FALL(I.K) 

15 CONTINUE 

CALL MATVEC(A.F.20,20.N,N,AF) 

CALL MATVEC(W0.AF.20O,2O.NS*N,N,W0AF) 
WRITE! 10. 156) !W0AF(I).I=1,8) 

156 format! 'W0AF«'/2(4E14. 4/)//) 

N0»0 

DO 20 I«1.NSN 

if!dabs(woaf!i)).lt.tol) N0=N0+1 
20 CONTINUE 

IF(NO. EO.NSN) GO TO 25 
ICAT(K)-1 
GO TO 60 
25 ICNT«0 

DO 30 I=1,N 
APAF!I)=AF(I) 

30 CONTINUE 
35 ICNT«ICNT+1 

CALL MATVEC!C.APAF. 10, 20. NS , N , CAPAF ) 
DO 40 I»1.NS 

CAAFN1!I.ICNT)«CAPAF!I) . 


40 


45 

50 


& 


:ONTINUE 

[F( ICNT. EO.N) GO TO 50 

:ALL MATVEC! AP , apaf . 20 . 20 . N , N , AP AF2 ) 

30 45 1*1. N 

APAF!I)*APAF2(I) 

:ONTINUE 
SO TO 35 

:all MAT1D!CAFNUL.20.N) 

:ALL 0RTRED(CAAFN1 ,NS.N. 10,20.CAFNUL,NS .o.cafnul. 

CAFSUB . IRKCAF . IC0DE3 . ZKNU 1 , WL . 10 . WF ) 
if! IRKCAF.NE. 1 ) GO TO 55 


FDF04410 

F0F04420 

F0F04430 

FDF04440 

F0FO4450 

F0FO446O 

F0FO4470 

FDF04480 

FDF04490 

FDF04500 

F0F04510 

FDF04520 

F0F04530 

F0F04540 

FDF04550 

FDF04560 

FDF04570 

FDF04580 

FDF04590 

F0F04600 

FDF046 10 

FDF04620 

FDF04630 

FDF04640 

FDF04650 

FDF04660 

FDF04670 

FDF04680 

FDF04690 

FDF04700 

FDF04710 

FDF04720 

FDF04730 

FDF04740 

FDF04750 

FDF04760 

FDF04770 

FDF04780 

FDF04790 

FDF04800 

FDF04810 

FDF04820 

FDF04830 

FDF04840 

F0FO485O 

F0FO486O 

F0FO487O 

F0F04880 

FDF04890 

FDF04900 

FDF04910 

FDF04920 

FDF04930 

FDF04940 

FDF04950 


/ 
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ICAT(K)«2 
GO TO 60 
55 ICAT(K)«3 
60 CONTINUE 

WRITE( 10. 166) ( ICAT(K) .K=1 .4) 

166 FORMAT( ' ICAT-' .413/) 

DO 65 K»1 .NF 

IF( ICAT(K).E0.3) RETURN 
65 CONTINUE 

CALL VMULFF(A.FALL.N.N.NF.20.20.AFALL.20.IER) 
WRITE( 10. 177) ( (AFALLd .«J).«J“1 .4) . I»1 .6) 
177 F0RMAT( ' AFALL- ' /6( 4E 14 . 4/ )// ) 

DO 68 1=1. NS 
DO 68 J=1.N 

CPF2(I.J)=CPF(I.d) 

68 CONTINUE 

DO 69 I«1.N 
DO 69 d-I.N 

AP2(I.d)»AP(I.d) 

69 CONTINUE 
IDIR=1 

C» IF AFALL2 IS INPUTTED TO PROGRAM. 

C= DELETE THESE NEXT 5 STATEMENTS 
NF2=NF 
DO 70 I=1,N 
DO 70 d=1 ,NF2 

AFALL2( I .d)-AFALL( I , J) 


70 CONTINUE ^ 

CALL SEPDET( AP2 .CPF2 . AFALL2 , N, NF2 , NS , I SEP . IDET .IDIR.C.CF.G. 

5 NU.NUM0D.ZKF.CPF.2MF) 

IF( IDET. EO. 1 ) RETURN 

CALL DGAIN(CPF2. AP2.EV,G.CF.CPF.NU,NUM0D.N.NF2.NS.0.0.0PR. 

6 ZKF.ZMF.AP.DPR) 

WRITE( 10. 127) ((DPR(I.d).J=1.4).I«1.6) 

127 F0RMATdDPR»'/6(4E14.4/)//) 

IDIR=0 

CALL SEPDET(A,C.FALL.N.NF.NS.ISEP,IDET.IDIR.C.CF, G.NU . NUMOO . 

5 2KF.CPF.ZMF) 

CALL DGAIN(C.A.EV.G.CF.CPF.NU,NUMOO.N.NF,NS.O, 1 .DPR. 

6 2KP;ZMF.AP.D) 


RETURN 

END 


• ORTRED PERFORMS ORTHOGONAL REDUCTION ON MATRIX V 


SUBROUTINE ORTRED( V , M1V. N 1 V . MV , NV , OMI ,NS , IC0DE4 , OM , 0MSU8 . 

& IRANK.IC0DE3.ZKNU1 .WL.IO.WF) 

INTEGER lO(IO).lCODEKlO) 

^DOUBLE PRECISIO . OMI (N1 V , N1 V ) . 0M( N1 V , N 1 V ) ^0MSUB(N1 V . N1 V ) . 

& VK(20).W(20).WL(20).WF(20. 10), 


, WTVK , ZKNU 1 . EPS . TOL . VKEPS 


DATA 10. EPS/8. 1 .0-4/ 

ICNT=0 

IRANK=0 

IC0DE2=0 


FDF04960 

FDF04970 

F0F04980 

FOF04990 

FDF05000 

F0F05010 

F0FO5020 

FDF05030 

F0F05040 

FDF05050 

FDF05060 

FDF05070 

F0F05080 

FDF05090 

FDF05100 

FDF051 10 

FDF05120 

FDF05130 

FDF05140 

FDF05150 

FDF05160 

FDF05170 

FDF05180 

FDF05190 

FDF05200 

FDF05210 

FDF05220 

FDF05230 

FDF05240 

F0F05250 

FDF05260 

FDF05270 

FDF05280 

FDF05290 

FDF05300 

FDF05310 

FDF05320 

FDF05330 

FDF05340 

FDF05350 

FDF05360 

F0F05370 

FDF0S380 

FDF05390 

FDF05400 

FDF05410 

FDF05420 

F0F05430 

FDF05440 

F0F05450 

FDF05460 

F0F05470 

FDF05480 

FDF05490 

FDF05500 
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205 



IC00E3-0 

FDF05510 


WRITE(I0,205) 

FDF05520 


F0RMAT( 'ENTERING ORTRED') 

FDF05530 


DO 5 I»1.NS 

FDF05540 


I0( I )*o 

FDF05550 


IC00E1(I)»0 

FDF05560 

5 

CONTINUE 

FDF05570 


DO 10 I»1 .MV 

FOF05580 


DO 10 U*1.NV 

FDF05590 


OM(I , J)*OMI(I , J) 

FDF05600 

10 

CONTINUE 

FDF05610 


DO 95 K=1 ,MV 

FDF05620 


T0L*0. 

FDF05630 


DO 20 0*1, NV 

FDF05640 


VK(0)*V(K.0) 

FDF05650 


VKEPS*VK(0)*EPS 

FDF05660 


T0L=DABS( DMAX 1( TOL . VKEPS ) ) 

FDF05670 

20 

CONTINUE 

FDF05680 


IF(TOL.LT.EPS) TOL*EPS 

FDF05690 


IC0DE2*IC0DE2+1 

FDF05700 


IARG=IC0DE2-( ( IC0DE2- 1 )/NS)*NS 

FDF05710 


IF(ICOOEKIARG).NE.O) go to 95 

FDF05720 


CALL MATVEC(0M, VK.NIV.NIV.NV.NV.W) 

FDF05730 


NZER0*0 

FDF05740 


DO 25 I=1.NV 

FDF05750 


IF(DABS(W( I ) ) . LT.TOL) N2ER0*NZER0+ 1 

FDF05760 

25 

CONTINUE 

FOF05770 


WTVK=0. 

FDF05780 


DO 29 1*1 .NV 

FDF05790 


WTVK*WTVK+W( I )*VK( I ) 

FDF05800 

29 

CONTINUE 

FDF05810 


IF(NZERO.NE.NV) GO TO 28 

FDF05820 


IC0DE1(IARG)*1 

FDF05830 


GO TO 95 

FDF05840 

28 

IF(DABS(WTVK) .GT.TOL) GO TO 30 

FDF05850 


IC00E1(IARG)»1 

FDF05860 


GO TO 95 

FDF05870 

30 

IF(ICODE3.NE.O) GO TO 35 

FDF05880 


1C0DE3=K 

FDF05890 

35 

DO- 40 I « 1 . NV 

FDF05900 


WL( I )*W(I) 

FDF05910 

40 

CONTINUE 

FDF05920 


IF( IC0DE4.NE. 1 ) GO TO 80 

FDF05930 


I0( IARG)*I0(IARG)+1 

F0F05940 


DO 60 I - 1 . NV 

FDF059S0 


WF(I,IARG)*W(I) 

FDF05960 

60 

CONTINUE 

FDF05970 

80 

IRANK-IRANK+1 

FDF05980 


WTVK«0. 

FDF05990 


DO 85 I»1,NV 

FDF06000 


WTVK=WTVK+W( I )*VK( I ) 

FDF06010 

85 

CONTINUE 

FDF06020 


ZKNU1*WTVK 

FDF06030 


DO 90 I»1 ,NV 

FDF06040 


DO 90 J*1 ,NV 

FDF06050 
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OMSUBd . J)*W( I)*W(0)/WTVK 
90 CONTINUE 

CALL MAT SUB ( OM , OMSUB , N 1 V , N 1 V , NV . NV . OM ) 

95 CONTINUE 
RETURN 
END 
C 

C* MATGEN GENERATES LARGER MATRIX C FROM A AND B 

SUBROUT'I NE MATGEN( A ,B,M1A,N1A,M1C,MA,NA.N, CSUB , CSU62 , C ) 
DOUBLE PRECISION A( M1 A . N1 A ) , B(N1 A , N1 A ) . C( M1C , N1 A ) . 

& CSUB(M1A,N1A),CSUB2(M1A,N1A),EPS 

DATA EPS/1. D-5/ 

DO 5 1*1 ,MA 
DO 5 J»1 ,NA 

IF(OABS(A(I, J)) .LT.EPS) A(I, J)*O.DO 
C(I . J)*A(I, J) 

CSUB(I,J)*C(I.J) 

5 CONTINUE 

IF(N.EO.I) RETURN 

NM1*N-1 

DO 15 K«1 ,NM1 

CALL VMULFF(CSUB ,B,MA,NA .NA,M1A,N1A .CSUB2.M1A, lER) 

DO 10 1-1 ,MA 
DO 10 0*1 ,NA 

CSUB(I.0)*CSUB2(I,J) 

C(MA*K+I.J)=CSUB(I,0) 

10 CONTINUE 
15 CONTINUE 
RETURN 
END 
C 

C* EVAS DETERMINES COEFFICIENTS P IN CHARACTERISTIC EQUATION 
C 

SUBROUTINE EVAS(N1 , N, EV , P ) 

REAL PDUM( 10) 

DOUBLE PRECISION P(N1),EV 
EVSNGL-SNGL(EV) 

DO 5 K*1 ,N 

BETA=( FLOArT( IFACT(N) )/( FLOAT( I FACT(N- I ) ) •FLOAT(IFACT( I ) ) ) ) 
& »(EVSNGL»*I) 

PDUM(K)»((-1.0)*«I)»BETA 

P(K)-08LE(PDUM(K)) 

5 CONTINUE 
RETURN 
END 
C 

C* MULTIPLIES MATRIX BY VECTOR 
C 

SUBROUTINE MATVEC(ARRAY,V.M1 ,N1 .M.N.RET) 

DOUBLE PRECISION ARRAY(M1 . N1 ) . V(N 1 ) . RET( M 1 ) 

DO 10 I-1.M 
RET(I )«0.D0 
DO 10 J»1.N 


FDF06060 

FDF06070 

FDF06080 

FDF06090 

FDF06100 

FDF061 10 

FDF06120 

F0FO6130 

F0F06140 

F0F06150 

FDF06160 

FDF06170 

FDF06180 

FDF06190 

FDF06200 

FDF06210 

FDF06220 

FDF06230 

FDF06240 

FDF06250 

FDF06260 

FDF06270 

FDF06280 

F0F06290 

FDF06300 

FDF06310 

FDF06320 

F0F06330 

FDF06340 

FOF06350 

FDF06360 

FDF06370 

FDF06380 

FDF06390 

FDF06400 

FDF06410 

FOF06420 

FDF06430 

FDF06440 

FDF0e4S0 

FDF06460 

FDF06470 

FDF06480 

FDF06490 

FDF06500 

FDF06510 

FDF0S520 

FDF06530 

FDF06540 

FDF06550 

FDF06560 

FDF06570 

F0F06580 

FDF06590 

FDF06600 
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RET(I)-RET(I)+ARRAY(I.J)*V(U) 

10 CONTINUE 
RETURN 
ENO 

■ AOOS MATRIX B TO MATRIX A 

SUBROUTINE MATAOO( A , B .Ml , N1 . M. N. RET ) 

DOUBLE PRECISION A(M1 ,N1 ) ,B(M1 ,N1 ) ,RET(M1 ,N1 ) 
00 10 I»1,M 
DO 10 J*1,N 

RET(I.J)»A(I.J)+B(I,J) 

10 CONTINUE 
RETURN 
ENO 

• SUBTRACTS MATRIX B FROM A 

SUBROUTINE MATSUB( A.B.MI ,N1 .M.N.RET) 

DOUBLE PRECISION A(M1 ,N1 ) ,B(M1 ,N1 ) , RET(M1 ,N1 ) 
DO 10 I»1,M 
DO 10 J=1 .N 

RET(I,d)-A(I.J)-B(I.J) 

10 CONTINUE 
RETURN 
END 

• GENERATES THE IDENTITY MATRIX 

SUBROUTINE MATID(A.N1 ,N) 

DOUBLE PRECISION A(N1.N1) 

DO 20 I « 1 . N 
DO 10 U»1 .N 
A(I.J)«0.00 
10 CONTINUE 

A(I.I)«1.0D0 
20 CONTINUE 
RETURN 
ENO 

• RAISES MATRIX Ar^fo POWER K 

SUBROUTINE MATP0W(A,N1 .N.K.AK2.AK) 

DOUBLE PRECISION A(N1 ,N1 ) . AK(N1 ,N1 ) , AK2(N1 ,N1 ) 
CALL MATID(AK.N1 ,N) 

IF(K.EO.O) RETURN 
00 10 IDUM-1.K 

CALL VMULFF(A.AK.N.N.N.N1 .N1.AK2.N1 . lER) 

DO 5 I>1 ,N 
DO 5 d“1.N 

AK(I.d)»AK2(l.d) 

5 CONTINUE 
10 CONTINUE 
RETURN 
ENO 


FDF06610 

FDF06620 

FOF06630 

FDF06640 

FDF06650 

FDF06660 

FDF06670 

FDF06680 

F0F06690 

FDF06700 

FDF06710 

F0F06720 

F0F06730 

FDF06740 

FDF06750 

FDF06760 

FDF06770 

FDF06780 

FDF06790 

FDF06800 

F0F06810 

FDF06820 

FDF06830 

FDF06840 

FDF06850 

FDF06860 

FDF06870 

FDF06880 

FDF06890 

FDF06900 

FDF06910 

FDF06920 

FDF06930 

FDF06940 

FDF06950 

FDF06960 

FDF06970 

FDF06980 

FDF06990 

FDF07000 

FDF07010 

FDF07020 

FDF07030 

F0F07040 

FDF07050 

FDF07060 

FDF07070 

FDF07080 

FDF07090 

FDF07100 

F0F071 10 

FDF07120 

F0F07130 

FDF07140 

FDF07150 
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• FINDS FACTORIAL OF N 

INTEGER FUNCTION IFACT(N) 
IFACT*1 

IF(N.LE. 1) RETURN 
DO 10 «J«2,N 
IFACT«IFACT*vJ 
10 CONTINUE 
RETURN 
END 


FDF07160 

FDF07170 

F0F07180 

F0F07190 

FDF07200 

FDF07210 

FDF07220 

FDF07230 

FDF07240 

FDF07250 

FDF07260 


+ + + 4 *+ 
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E 3 X 3 3 S 3 
: 3 3 3 X 3 3 


K 3 3 3 3 3 S 


I PC 
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10 3 3 

- 15.0E0 
1.0 
O. 

O. 

O. 

O. 

O. 

O. 

O. 

O. 

O. 

o. 

o. 

o. 

0 . 

0 . 

o. 

0 . 

1 .0 

o. 

rS2.6 O. 

0 . 

O. 

O. 

O. 

0 . 

0, 

O. 

O. 

o. 

o. 

0 . 

o. . 


00 o. 

04 O. 

OO o. 

.30 O. 

00 o. 

41 O. 


0 1 

o. 0 . 

0 . 0 . 

O. 0. 

0, O. 

O. 1-0 

0 , o. 

o. o, 

o. o. 

0 . 0 . 

o. o. 

o. o. 

0 . o. 

o. 0 . 

0 . o. 

0 . o. 

o. 0 . 

0 . 0 . 

o. o. 

O. 0 . 

o. o. 

0 . 0 . 

O. 1.0 

o. o. 

-10181. O. 

o. o. 

0 . 0 . 

o. o. 

o. 0 . 

o. o. 

o. 0 . 

0 . 0 . 

O. 0 . 

-1.59 0. . 

-0,753 O. 

- 1.01 o. 

1.09 o. 

0.00 o. 

0.00 o. 


o. o. 

0 . 0 . 

0 . o. 

0 . 0 . 

0 . o. 

o. o. 

0. O. 

0 . 0 . 

O. 1-0 

0 . o. 

-130.38 O. 

o. o. 

0 . o. 

0 . 0 . 

o. o. 

o. o. 

0 . 0 . 

o. o. 

0 . o. 

0 . o. 

0 . o. 

o. O. 

o. o. 

o. o. 

o. o. 

0 . 1.0 

o. o. 

-22556. O. 

o. o. 

o. o. 

0 . 0 . 

0 , o. 

-1.60 O. 

-0.465 O, 

-0.123 0. 

0.356 O. 

1.21 O. 

1.40 0. 


O. O. 

0 . 0 . 

0. O. 

O. 0. 

0 . 0 . 

O. 0. 

0 . 0 . 

0 . o. 

O. 0. 

o. 0 . 

0. O. 

0 . 0 . 

0 . 1.0 

0 . o. 

-983.45 0. 

0 . 0 . 

0 . 0 . 

0 . 0 . 

O. 0. 

0 . 0 . 

O. 0. 

0 . 0 . 

0 . 0 . 

O. 0. 

0 . 0 . 

O. 0. 

0 . 0 . 

O. 0. 

0. O. 

o. 1-0 

o. 


-43683. 
1.321 
-O. 181 
-0.876 
-0.553 
0-00 
0.00 


O. O. 

1.992 1.992 1.992 

O. 0. O. 

-3.167 -2.012 0.000 

0. O. O. 

-3.187 -0.2450 2.410 

0. O. O. 

2.631 - 1 . 745 0.000 

0. 0. O. 

2.072 -2.590 2.809 

0. O. 0. 

-1.500 2.171 0.000 

O. O. O. 

-0.9263 0.7092 2.789 

O. 0. O. 

-0.3606 -1.102 0.000 
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APPENDIX B 


FILE: FDSIM FORTRAN A 


VM/SP CONVERSATIONAL MONITOR SYSTEM 


C* 

C 

C 

C 

C 


PROGRAM FDSIM SIMULATES CONTINUOUS 
SYSTEM AND FILTER DYNAMICS 
FOR FAILURE DETECTION AND IDENTIFICATION 








C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


INPUT: NSYS - NUMBER OF SYSTEM STATES 

NFIL - NUMBER OF FILTER STATES 
NS - NUMBER OF SENSORS 
NA - NUMBER OF ACTUATORS 
NF - NUMBER OF FAILURE EVENTS 
IAS - («1) for actuator FAILURES 
(«2) FOR SENSOR FAILURES 
OMU - ACTUATOR INPUT FREQUENCY 
ISR - SAMPLING RATE (HZ) 

NTSEC - INTEGRATION FREQUENCY (HZ) 

IFAIL • ^ Of FAILED ACTUATOR OR SENSOR 
TFAIL - TIME OF FAILURE (SEC) 

TLAST - TIME OF SIMULATION 
A - SYSTEM MATRIX 

BS - CONTROL EFFECTIVENESS MATRIX FOR SYSTEM 

CS - MEASUREMENT MATRIX FOR SYSTEM 

CF1 - MATRIX OF OUTPUT ERROR DIRECTIONS FOR 

ACTUATOR FAILURES (COLUMNS ARE VECTORS) 
D - DETECTOR GAIN 
X - INITIAL STATES 

OUTPUT: RE - TRANSFORMED OUTPUT ERROR VECTOR 








DIMENSION A(20.20).BS(20. 10).BF(20, 10).CS( 10.20).U( 10). 

& 0(20. 10).CF1( 10. 10).X(40) 

& R( 10. 10) .CFTCF( 10. 10) .CFTCF ( 10. 10 .CFTTK 10. 10) . 

& WK(40).WK1(40.9).RE(10).C1A(24).CF{ 10.20). 

& XS(20).XF(20).YS(10).YF(10), 

8. TIME(5000).RE1(5000).RE2(5000).RE3(5000) 

EXTERNAL FCN . ^ 4 c a a^/ 

DATA IN.IO.INO.TOL.EPS.NW/9.10.1.0.001.1.E-4 40/ 

COMMON A.BS.BF.CS.CF.D.E.U.OMU.NSYS.NFIL.NA.NS 

READ( IN?5)^ NSYS ! NFIL !ns!nA.NF. IAS. OMU, ISR. NTSEC. IFAIL. TFAIL. 

& ((A(I.J).J«1.16).I-1.16), 

5 ((BS(I.d).J»1.8),I»1.16). 

a ((CS(I.d).0-1.16).l-1.S). 

6 ((CF1(I,d).d«1.3).I“1.8), 

& ((D(I.J).J"1.3).I«1.16). 

5^F0RMAT(6I3/F1o!4/I3/213!2F1O,4/32(8F10.y).16(8F1^ 

& 16(8F10.4/).8(3E16.6/).16(3E16.6/).4(8F10.4/)) 

WRITEdO.7) NSYS.NFIL.NS.NA.NF .NTSEC , I FAIL , TF A I L. TLAST . 
a ((CFl(I.d).d*1.4).I-1.4), 


FDS00010 
FDS00020 
FDS00030 
FDS00040 
FOS00050 
FDS00060 
F0S00070 
FDS00080 
FOS00090 
F0SOO1OO 
FDS001 10 
FDS00120 
FDS00130 
FDSO0140 
FDS00150 
FDSO016O 
FDS00170 
FDS0O180 
FDS00190 
FDS00200 
FDS00210 
FDS00220 
FDS00230 
F0S00240 
FDS00250 
FD500260 
FDS00270 
FDS00280 
FDS00290 
FDS00300 
FDS00310 
♦FD500320 
FDS00330 
FD500340 
FDS00350 
FDS00360 
F0S00370 
FDS00380 
F0S0O39O 
FOS00400 
FOS00410 
FDS00420 
FDS00430 
F0SOO44O 
TLAST. FDS00450 
F0S00460 
FDS00470 
F0S00480 
FDS00490 
F0S00500 
FDS00510 
FDS00520 
F0S00530 
FOSOO540 
FOS00550 
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FILE: FDSIM 


FORTRAN A 


VM/SP CONVERSATIONAL MONITOR SYSTEM 


7*F0RMAT('NSYS-'?12!''NFIL«^I2.5X 'NS-'.I2.5X. 'NA-'^I2^5X^^^ 

& 'NTSEC-' .I3.5X. 'IFAIL*' .I2.5X. 'TFAIL- .F8.3.5X. TLA5T . 

& F8. 3/'CF1«'/4(4E14.4/)/'D«'/4(4E14 .4/)/) 

NE«NSYS+NFIL 
DO 10 I*1,NFIL 
DO 10 J”1,NA 

8F(I.d)«8S(I.d) 

10 CONTINUE 

DO 12 1*1. NS 
DO 12 J»1,NFIL 
CF(I , J)-CS(I.J) 

12 CONTINUE 

DO 20 I * 1 , NS 
00 15 J»1.NF 
S(I .vJ)"0. 

15 CONTINUE 

S(I,I)“1.0 

20 CONTINUE 
IF(IAS.E0.2) GO TO 21 

CALL VMULFM(CF1 .CF1 .NS.NF.NF. 10. 10.CFTCF . 10. lER) 

CALL LINV1F(CFTCF.NF.10.CFTCFI.0.WK lER) 

CALL VMULFP(CFTCFI .CF1 .NF .NF .NS. 10. 10.CFTTI , 10. lER) 

CALL VMULFF(S.CFTTI.NS.NF.NS. 10. 10.R. 10. lER) 

21 C0SGAM=CF1(2.1)/S0RT{CF1(2.1)**2+CF1(3,1)**2) 

5INGAM=CF1(3. 1 )/SORT(CF 1 (2 . 1 )»*2+CF 1 ( 3 . 1 ) * *2 ) 

R( 1 , 1 )*1 .0 
R( 1 .2)*0. 

R( 1 .3)»0. 

R(2. 1 )»0. 

R(2.2)=C0SGAM 
R(2.3)*SINGAM 
R(3. 1 )*0. 

R(3.2)«-SINGAM 

R(3.3)-C0SGAM 

24 T-0. 

CALL FCN(NE,T.X.XD0T) 

WRITE(I0,25) , 

25 F0RMAT(2X.'TIME'.16X.'E1'.12X.'E2'.12X.'E3'/) 

WRITE(I0.3O) T.(RE(I).I*1.3) 

30 F0RMAT(2X.F6.3. 10X.3E14.4) 

ILAST*INT(TLAST)*NTSEC 

IPLOT*ILAST 

NFE*NTSEC/ISR 

IN0EXE*0 

IWRT-0 

00 50 K»1 . I LAST 

TEND»FLOAT(K)/FLOAT(NTSEC) 

C* NEXT STATEMENTS FOR DATA-SAMPLING ONLY 
C INDEXE-lNDEXE+1 

C IF(INDEXE.NE.NFE) GO TO 38 


C NP1»NSYS+1 

C DO 32 I»1.NSYS 

C XS(I)-X(I) 


FDS00560 

FDS00570 

FDS00580 

F0S00590 

FDS00600 

FDS00610 

F0S0O620 

FDS00630 

FDS00640 

FDS00650 

FDS00660 

FDS00670 

FD500680 

FDS00690 

FDS00700 

FDS007 10 

FDS0O720 

FDS00730 

FDSOO740 

FDS00750 

FDS00760 

FDS00770 

FDS00780 

FDS00790 

FDS00800 

FDS00810 

FDS00820 

FDS00830 

FDS00840 

FDS00850 

FDS00860 

FDS00870 

FDS0C880 

FDS00890 

FDS00900 

FDS00910 

FDS00920 

FDS00930 

FDS00940 

FDS00950 

FDS00960 

FDS00970 

FDS00980 

FDS00990 

FDS01000 

FDS01010 

FDS01020 

FDS01030 

FDS01040 

FDS01050 

F0S01060 

FDS01070 

FDS01080 

F0S01090 

FDS01 100 
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FILE: FDSIM 


fortran a 


VM/SP CONVERSATIONAL MONITOR SYSTEM 


32 


34 


38 


40 

415 

416 

41 


45 


CONTINUE 
00 34 I'NPl.NE 
XF(I-NSYS)»X(I) 


CALL MATVEC(CS.XS. 10 , 20.NS . NSYS . YS 

CALL MATVEC(CF,XF. 10.20.NS.NFIL.YF 
CALL VECS(YS.YF. 10.NS.E) 

INOEXE-0 

IF(T.LT.TFAIL) GO TO 41 
IF(IFLAG.EO. 1) GO TO 41 
IFLAG*1 

IF(IAS.E0.2) GO TO 415 
DO 40 I»1,NSYS 
BS(I . IFAIL)«0. 


) 

) 


CONTINUE 


GO TO 41 
DO 416 J«1.NSYS 
CS(IFAIL.J)«0. 


CALL DVERK(NE.FCN.T.X.TEND.T0L.IND,C1A.NW,WK1 ,IER) 

IF(IND.LT.O .OR. lER.GT.O) GO TO 62 

CALL MATVEC(R.E. 10. 10.NS.NS.RE) 

TIME(K)»T 
RE1 (K)«RE{ 1 ) 

RE2(K)-RE(2) 

RE3(t<)»RE(3). 

lWRT-IWRT+1 

IF( IWRT.NE.5) GO TO 50 
WRITE(I0,45) T. (RE(I) . 1=1 .3) 

F0RMAT(2X,F6.3. 10X.3E14.4) 


IWRT»0 
50 CONTINUE 

CALL PICTUR(a.0,3.0. 'TIME(SEC) 
a TIME.REI . IPLOT.O. , 

a TIME.RE2. IPLOT.O. . 

a TIME. RE3, IPLOT.O. . 

CALL ENDPLT( 12.0.0.0.999) 

GO TO 100 
62 WRITE(I0.75) 

75 FORMAT( ' IND<0 OR IER>0') 

100 STOP 
END 


' .9. ' 

O, 


1 . 

2 ) 


error(ft) 


.9. 


C 

C 

C 


IIMENSION J<NE) 

E( 10), 0(20, 10), A( 20, 20) , AXS(20) , AXF ( 20) . BSU( 20) . 
BFU(20) ,0E(20) , AXFBFU(20) .CF( 10,20) 

:OMMON a.bs.bf.cs.cf.d.e.u.omu.nsys.nfil.na.ns 
J( 1)«EXP(-0.2»T)»SIN(0MU»T) 

J(2)-EXP(-O.2*T)»SIN(0MU»T) 

3)«EXP( -O. 2*T)»SIN(0MU«T ) 


NP1-NSYS+1 


F0S01 1 10 

FDS01 120 

FDS01 130 

FOS01 140 

F0SO1150 

F0S01160 

FDS01170 

F0S01 180 

FDS01190 

FDS01200 

F0S01210 

F0S01220 

F0S01230 

FDS01240 

FDS01250 

FDS01260 

FDS01270 

FDS01280 

FDS01290 

FDS01300 

FDS01310 

FDS01320 

FDS01330 

FDS01340 

FDS01350 

FDS01360 

FDS01370 

F0S01380 

FDS01390 

FDS01400 

FDS01410 

FDS01420 

FDS01430 

FDS01440 

FDS01450 

FDS01460 

FDS01470 

FDS01480 

FDS01490 

FDS01500 

FDS01510 

FDS01520 

FDS01530 

F0S01540 

FDS01550 

FDS01560 

FDS01570 

FDS01580 

FDS01S90 

F0S01600 

FOS01610 

F0S01620 

F0S01630 

FDS01640 

F0S01650 
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FILE: FDSIM 


FORTRAN A 


VM/SP CONVERSATIONAL MONITOR SYSTEM 


DO 5 I-1.NSYS 
XS(I)-X(I) 

XSD0T(I)-XD0T(1) 

S CONTINUE 

DO 10 I«NP1,NE 
XF(I-NSYS)«X(I) 

XFDOT( I -NSYS)-XDOT( I ) 

10 CONTINUE 

CALL MATVEC( A . XS . 20. 20, NSYS . NSYS . AXS ) 
CALL MATVEC(8S.U,20. 10, NSYS ,NA ,8SU ) 
CALL VECP(AXS.8SU.20.NSYS.XSD0T) 

CALL MATVEC(A,XF.20.20.NFIL.NFIL.AXF) 
CALL MATVEC(CS.XS. 10. 20. NS , NSYS . YS ) 
CALL MATVEC(CF,XF , 10. 20. NS . NFIL , YF ) 
CALL VECS(YS. YF. 10.NS.E) 

CALL MATVEC(0,E.20. 10. NFIL. NS. DE) 
CALL MATVEC(8F .U. 20, 10, NFIL . NA . BFU ) 
CALL VECP( AXF,8FU,20,NFIL. AXFBFU) 
CALL VECP( AXFBFU.de. 20, NFIL, XFDOT) 

DO 15 1*1. NSYS 
XDOT(I)=XSDOT(I) 

15 CONTINUE 

DO 20 I=NP1,NE 

XDOT( I )=XFDOT( I-NSYS) 

20 CONTINUE 
RETURN 
END 


SUBROUTINE MATVEC( ARRAY . V. Ml .NI.M.N.RET) 
DIMENSION ARRAY(M1 ,N1 ) ,V(N1 ) .RET (Ml ) 

DO 10 1=1, M 
RET( I )»0. 

DO 10 J»1.N 

RET( I )»RET( I )+ARRAY( I . J)*V(d) 

10 CONTINUE 
RETURN 
END 


SUBROUTINE VECP( VI , V2 . N1 . N. VADO ) 
DIMENSION V1 (N1 ) . V2(N1 ) . VAD0(N1 ) 
DO 5 1*1, N 

VAD0(I)«V1(I)+V2(I) 

5 CONTINUE 
RETURN 
END 


SUBROUTINE VECS( VI , V2 . N1 , N, VSU8 ) 
DIMENSION V1(N1).V2(N1).VSUB(N1) 
DO 5 I»1,N 

VSUB(I)»V1(I)-V2(1) 

5 CONTINUE 
RETURN 
END 


FDS01660 

F0SO1670 

FDS01680 

F0S01690 

F0S01700 

F0S01710 

FDS01720 

FDS01730 

F0S01740 

FDS01750 

FDS01760 

FDS01770 

F0S01780 

FDS01790 

FDS01800 

FDS01810 

FDS01820 

FDS01830 

FD501840 

FDSO1850 

FDS01860 

FDS01870 

FDS01880 

FDS01890 

FDS01900 

FDS01910 

FDS01920 

FDS01930 

FDS01940 

FDS01950 

FDS01960 

FDSO1970 

FDS01980 

FDS01990 

FDS02000 

F0S02010 

FDS02020 

FDS02030 

FDS02040 

FDS02050 

FDS02060 

FDS02070 

FDS02080 

F0SO2O9O 

F0S02100 

F0S021 10 

FDS02120 

F0S02130 

F0S02140 

F0SO215O 

FDS02160 

F0S02170 

F0S02180 

F0S02190 

FDS02200 

FDS02210 

F0S02220 

FDS02230 

FDS02240 
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FILE; FDSIMS DATA 


VM/SP CONVERSATIONAL MONITOR SYSTEM 


16 10 
20.0 
64 


3 3 1 


192 2 

1.0 2 

.01 






0. 

1.0 

0. 

0. 

0. 

0, 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0, 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

1.0 

0. 

0. 

0. 

0. 

0. 

O. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0, 

0. 

0, 

0. 

0. 

0. 

0, 

0. 

0. 

0. 

0. 

0. 

1.0 

0. 

0. 

0. 

0, 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

-130.38 

0. 

0. 

0. 

0. 

o. 

0. 

O. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

1.0 

0. 

0. 

0. 

0. 

0. 

0, 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

-983.45 

0. 

0. 

o. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

1 .0 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

-3752.6 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

1 .0 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

O. 

0. 

-10181 . 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

1 .o 

o. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

-22556. 

0. 

0, 

0. 

0. 

0. 

0. 

0. 

0, 

0. 

0. 

0. 

0. 

0. 

0, 

0. 

0. 

0. 

0. 

1 .0 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

-43683. 

0. 

0. 

0. 

0. 






1.992 

1.992 

1.992 






0. 

0. 

0. 






-3. 167 

-2.012 

0,000 






0. 

0, 

0. 






-3. 187 

-0.2450 

2.410 






0. 

0. 

0. 






2.631 

-1.745 

0.000 






0. 

0. 

0. 






2.072 

-2.590 

2.809 






0. 

0. 

0. 






-1.500 

2.171 

0.000 






0. 

0. 

0. 






-0.9263 

0.7092 

2.789 






O. 

0. 

O. 






-0.3606 

-1.102 

0.000 




1.321 


1 .00 

0. 

-1.59 

0. 

-1.60 

0. 

0. 

1 .04 

O. 

-0.753 

0. 

-0,465 

0. 

-0. 181 

0. 

1 .00 

0. 

-1.01 

0. 

-O. 123 

0. 

-0.876 

0. 
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FILE: 

FDSIMS 

DATA A 


-1.30 

0. 

1.09 

0. 

1 .00 

0. 

0.00 

0. 

1 .41 

0. 

0.00 

0. 


VM/SP CONVERSATIONAL MONITOR SYSTEM 


0.356 

0. 

-0.553 

0 

1,21 

0. 

0.00 

0 

1.40 

0. 

0,00 

0 


0, 177572E+02 
0.584315E+00 
O. 105725E+01 


0.584335E+CXD 
0.894988E+01 
-O. 195635E+01 


O. 105736E+01 
-O. 195613E+01 
0.886879E+01 


0.214951E+01 
O. 193940E+02 
-0.560015E+01 
-0.383867E+02 
-0.586796E+01 
-0.207824E+02 
0.575753E+01 
-0.884495E+02 
0.294047E-01 
-0.484272E+03 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 

0. O. 

0 . 0 . 

0 . 0 . 

0 . 0 . 


0. 149971E+02 
0.635415E+02 
0.688610E+01 
-0.240117E+01 
0.746371E+01 
0.575569E+02 
-O. 165995E+02 
-0.434314E+O3 
0.562981E+02 
0.569748E+03 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 

0 . 0 . 0 . 

0 . 0 . 0 . 

0 . 0 . 0 . 

0. O. O. 


0.981505E+01 
0.627022E+02 
-0.487336E+01 
-0.486003E+02 
O. 1 19337E+01 
O. 1 13480E+02 
-0.985363E+01 
0.587169E+02 
0.362257E+01 
O. 1 17701E+04 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 


O, O. 

0 . 0 . 

0. O. 

0 . 0 . 
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